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ABSTRACT 


A general combustion and heat transfer model for porous 
media subject to Darcy flow is formulated. The transient 
one dimensional model treats the combustion process in two 
phases. During the initial phase, combustion occurs within 
the porous medium. The second phase occurs when the exo- 
thermic reaction moves to the air inlet surface of the medium 
resulting in surface recession. The temperature dependency 
of the system parameters and thermophysical properties is 
taken into account. An analysis of combustion in a carbon 
porous medium is presented, as well as an assessment of the 


accuracy of the model. 
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I. INTRODUCTION 


Studies were conducted at the Naval Weapons Center, China 
Lake, California to assess the survivability of the F-18 
aircraft when exposed to open pool fires. A portion of the 
F-18 aircraft is constructed of graphite-epoxy materials. 

When this composite material is exposed to a severe thermal 
environment, the epoxy burns off at a low temperature (about 
400-500 degrees Fahrenheit) leaving behind a porous graphite 
mat. Experiments were conducted by J. Fontenot [1] to evaluate 
the combustion behavior of the porous graphite or carbon mat. 
Due to the complex behavior of the combustion observed for the 
carbon fibers, a mathematical model to simulate this behavior 
was developed by Vatikiotis [2]. The model was formulated to 
consider a porous medium composed of cylindrically shaped 
fibers, typical of fiber reinforced composite materials. The 
details of the initial effort were presented by Vatikiotis. 
Significant improvements have been made to the initial combus- 
tion model, and the model has been extended to consider 
Spherically shaped particles, typical of thermal flow reactors. 
The initial work of Vatikiotis [2], where it applies to the 
present investigation is included in the discussion for the 
sake of completeness. Other analytical investigations taking 
a Similar approach, but concerned with other aspects of com- 
bustion or thermally active porous media, have been performed 


by Kordylewski [3], Sahota and Pagni [4], and Mehta, Sams and 
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Luss [5]. Kordylewski performed a numerical analysis of homo- 
geneous combustion in a two dimensional reactor. The objec- 
tive of Kordylewski's investigation was to show that the flow 
within the reactor influenced the conditions necessary for 
ignition. Sahota and Pagni considered the transient behavior 
of heat transfer in a porous medium when exposed to a fire. 
Their analysis was part of an overall investigation for deter- 
Mining the internal stresses in a structure caused by a fire. 
Combustion of the porous medium was not a consideration. 
Mehta, Sams, and Luss investigated the transient temperature 
response of a packed-bed reactor when the entering fluid 
temperature was decreased. Their model assumed that the 
temperatures of the porous solid and the fluid were the same. 

Slattery [6], Yaron [7], and Kassoy [8] present three 
approaches for mathematically modelling a porous media. 
Slattery's and Yaron's approaches were similar in using inte- 
gral methods in their development. Slattery considered only 
mass transfer whereas Yaron treated transport processes in 
general. Kassoy's approach was based on performing Taylor 
series expansions of energy transport processes with respect 
to a differential volume of a porous medium. An energy balance 
yielded two differential heat transfer equations, one for the 
porous solid and one for the fluid. A fundamental requirement 
of this approach is that the particle diameters remain small 
with respect to the thickness of the porous medium. 

Although the above investigations were concerned with 


Specific aspects of thermally active porous media, the present 


is 





combustion model was developed to treat a larger class of 
problems (1.e., combined mass transfer, heat transfer, and 
combustion). In contrast to other investigations, the present 
model takes into account the temperature dependency of all 

the thermophysical properties of the air and the geometric 
parameters of the porous medium. The physical parameters such 
as porosity, permeability, and thickness affect the magnitude 
and behavior of pressure-driven air flow within the porous 
medium. As discussed by Kordylewski [3], the air flow influ- 
ences the conditions for lgnition. Therefore, sensitivity 
analyses were performed using the combustion model showing 

the effects of these parameters on combustion. In addition, 

a functional relation between the internal pore velocity and 
an initial uniform temperature distribution needed to sustain 
combustion is presented. The initial work of Vatikiotis [2] 
showed that the initial conditions are important in deter- 
mining whether or not sustained combustion will occur. An 
analysis is presented showing the dependence of combustion 

on the shape of the initial temperature. In addition, the 
effects of the boundary conditions are examined using the 
combustion model. Specifically, there is a change in the 
behavior of combustion when insulated boundaries on the porous 
solid are changed to permit heat transfer by radiation. The 
Semenov model is used to explain the behavior observed during 
the analyses discussed above. Although the Semenov model 


was originally proposed for homogeneous combustion, the 
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results show that the Semenov model can be applied to hetero- 


geneous combustion in porous media. 
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Pie DeoekirLTLlON OF THE PROBES 


The problem under investigation is that of a porous slab, 
with a pressure differential across its thickness, in a severe 
thermal environment. The effect of the pressure differential 
is to induce air flow through the porous medium. This internal 
air flow produces two opposing effects, (1) internal convection 
heat transfer, and (2) a supply of oxygen for heat generation 
by combustion. Whether the porous medium moves towards sus- 
tained combustion or extinguishment depends on the interaction 
of these effects. A mathematical model was formulated to 
provide an understanding of this interaction and its effect 
on thermal behavior. 

The mathematical model was formulated as follows. Energy 
balances on the porous solid and the convected air provide 
heat transfer equations for each. The heat transfer mechan- 
isms included in the model are (1) conduction, (2) convection, 
and (3) radiation. In addition, nonvolatile combustion 1s 
included in the porous solid heat transfer equation as a heat 
generation term of Arrhenius type. 

The conservation of species law was applied to the oxygen 
molecule concentration. The resulting mass transfer equation 
includes the transfer mechanisms of (1) molecular diffusion, 
and (2) species transport by convection. A term accounting 
for the consumption of oxygen due to combustion is included 


as well. Darcy's law was taken as the constitutive equation 
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for flow through the porous medium. The pressure gradient, 
needed to calculate pore velocity, is obtained by solving a 
combined Darcy's law and continuity equation. 

As the behavior of the system occurs over a large tempera- 
ture range, all thermophysical properties which depend upon 
temperature, such as conductivity, viscosity, and air density, 
are treated as temperature dependent properties. Moreover, 
as combustion consumes the porous solid, changes in porosity 
and permeability are accounted for in the transient analysis. 

In order to provide a general model, a number of boundary 
conditions are included in the formulation of the model. A 
detailed discussion of boundary conditions is presented in 
Section III.H. The final system of equations to be solved are 
four tranSient, nonlinear, coupled partial differential equa- 
tions. The three heat and mass transfer equations were solved 
by a Galerkin formulation of the finite element method. The 
remaining combined Darcy's law-continuity equation was solved 
by a shooting method. The details of the solution procedure 


are presented in Appendix C. 
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IIIf. THEORY AND BACKGROUND 


A. DESCRIPTION OF THE POROUS MEDIUM 

For this investigation, a porous medium is defined as a 
solid containing interconnected pores (i.e., path for airflow). 
There are two classes of porous media, consolidated and uncon- 
solidated. Consolidated porous media are those where the 
solid constituent or matrix remains rigid. Consolidated por- 
ous media may be comprised of either individual particles 
attached to one another (e.g., sintered metal), or a solid 
where a portion has been removed (e.g., charcoal). Unconsoli- 
dated porous media are comprised of discrete particles (e.g., 
granular beds). The properties of an unconsolidated porous 
medium can change if the porous medium De agitated or settling 
takes place. Both consolidated and unconsolidated porous 
media can have spatially varying properties. 

The physical properties which characterize porous media 
(1.e., allow comparison of porous media without regard to class 
or form) are (1) porosity, (2) specific internal area, (3) pore 
Size or diameter, and (4) tortuosity. These properties are 
common to both consolidated and unconsolidated porous media. 
Porosity, p, is defined as the ratio of void volume to total 
volume. The specific internal area, z, is the ratio of in- 
ternal surface area to bulk volume. The dimension is the 
reciprocal of length. Pore diameter is used to describe the 


pore system of a porous medium. The actual shape and size of 
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a pore may be irregular and complex, and would be difficult 
to describe geometrically. In practice, there are many con- 
ventions used to define the pore diameter. For this inves- 
tigation, a hydraulic diameter theory is used and is discussed 
later in this section. The tortuosity of a porous medium is 
defined as the ratio of the length of the flow path of a fluid 
particle to the straight line distance. Originally, tortu- 
osity was considered a kinematic property, but was subsequently 
adopted for characterizing porous media. Tortuosity is also 
discussed later in this section. Scheidegger [9] discusses 
various methods used to measure properties of porous media. 
Though most methods discussed are based on experimental tech- 
niques, the above properties for this investigation are calcu- 
lated based on an idealization of the geometry of a porous 
medium. 

The porous medium was modelled as shown in Figure III.1l. 
In Figure III.l, D is the distance between the fiber or parti- 
cle centers, and dis the particle diameter. The geometry is 
idealized since the actual porous medium may be irregular in 
particle diameter and distribution. For cylindrical fibers, 
the porosity, which was defined as the ratio of void per unit 


volume is, 
TT 2 
p = 1- g(d/D) Geta. 1) 


and for spherical particles, 
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Idealized geometry of a porous medi 


FIGURE III.1 





pae= OL = 7(d/D) ° (Tome. ) 


The pore diameter is obtained by an expression proposed by 


@eeman [10], 
§ = 4p/z (III. 3) 


where z is the specific internal area or the particle surface 
area per unit volume of porous medium. Expression III.3 is 
analogous to the more familiar form of mean hydraulic diam- 
eter, 4V/A, where V is the void volume and A is the wetted 
surface area. The specific internal area, z, based on the 
idealized geometry of Figure III.1 is calculated by, 


ae = a/D* (IIT. 4) 


for cylindrical fibers, and by, 


Zz = =n d°/D (III.5) 
for spherical particles. Expressions III.4 and III.5 assume 
one half the total internal surface area is effective for 
convection heat transfer. This fractional amount of total 
area was an estimate based on Fontenot's experimental results 


and does not apply to porous media in general. A more general 


approach, as presented by Scheidegger [9], is the Kozeny 
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equation given by, 
z = (sp°ym) +77 (OiiaienG ) 


In equation III.6, m is the specific permeability or hydraulic 
conductivity (discussed in Section III.B), and s is the Kozeny 
constant. The Kozeny constant is usually taken as .2 for 
specific surface calculations. Advantages to using the Kozeny 
equation are (1) the calculated values of specific internal 
area are in fair agreement with the experimental values, and 
(2) the calculation is independent of particle shape. A 
disadvantage is that the Kozeny equation fails for highly 
porous fibrous media. The tortuosity, tT, is the ratio of the 
length of the flow path for a fluid particle to the straight 
line distance. For the geometric configuration of Figure III.1l, 
the tortuosity depends on the ratio d/D. Carman [10] pre- 
sents a table of measured tortuosity factors of various materials 
and geometries, and points out the differences between the 
analytical determinations of tortuosity. Here, we adopt his 
recommended value of 1.4. Since particle size decreases as 
the carbon is consumed, all geometric properties which depend 
on fiber or particle diameter are functions of time and posi- 
tion. The model assumes that the carbon matrix remains rigid 
as the particle diameter decreases, and thus, porosity in- 
creases with combustion. In his discussion on the packing 
theory of spheres, Scheidegger [9] states the highest reported 


porosity for spherical particles in a stable configuration is 
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.875. Carman [10] reports of investigations of fibrous porous 
media with porosities as high as .99. The change in particle 
diameter is accounted for, and will be discussed in the sec- 
tion on the Arrhenius expression (Section III.D) for carbon 


combustion. 


B. DARCY'S LAW AND PORE VELOCITY 


Reynolds number for porous media is defined by 


Rew = (9 ud)/u (aie?) 


where u is the local pore velocity, P. 1s the mass density of 
air, and yp is the dynamic viscosity. Depending on the magni- 
tude of the Reynolds number, the motion of the fluid may be 
dominated by molecular, viscous, or inertial effects. Most 
investigations of fluid flow in porous media have been in the 
range of Reynolds number where the flow is dominated by vis- 
cous and inertial effects. The Navier-Stokes equations are 
applicable in describing the fluid motion in this range of 
Reynolds number. However, because of the convoluted geometry 
and the necessary boundary conditions (i.e., u = 0 at a solid- 
fluid interface), solution of the Navier-Stokes equations are 
difficult for a porous medium. Extensive experimental work, 
as discussed by Scheidegger [9] has shown that fluid flow in 
porous media is governed by Darcy's law for the range of 
Reynolds number where viscous effects dominate. The upper 


limit (high velocity end) Reynolds number reported in the 
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experimental investigations varied from .1 to 75. As Reynolds 
number increases, inertial effects increase. However, Boffa 
[11] has shown that for a fixed Reynolds number, inertial 
effects diminish with increasing air temperature. For the 
problems covered in this investigation, the Reynolds number 
did not exceed a value of 5. 

Neglecting body forces, Darcy's law for one dimensional 


flow is, 


where Q is the filter velocity or the volumetric flow rate per 
unit of cross-sectional area, and dP/dx is the pressure gradi- 
ent. The specific permeability or the hydraulic conductivity 


of the porous medium is defined by, 
2 
m = 96p (6/tT) (III.9) 


Bepression III.9 is based on a capillaric-serial model given 
by Scheidegger [9]. It should be noted that, in a physical 
sense, permeability and porosity are not related. Porosity 
1s a quantifiable property of the porous medium. Whereas, 
permeability is a constitutive property (1.e., a property 
specified by a constitutive equation, Darcy's law). As in 
the case of Fourier's law of heat transfer, Darcy's law is 


not derived from first principles, but has been obtained 
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through exhaustive experimental analyses. The Dupuit- 
Forcheimer assumption, presented by Carman [10], relates 


the local pore velocity to the filter velocity by, 
Q = pu ee. 20) 


The hypothesis of the Dupuit-Forcheimer assumption is that the 
local pore velocity is greater than the filter velocity. 
Noting that the actual velocity in a single pore is not a 
constant, but rather a function of the location within the 
pore, the Dupuit-Forcheimer assumption defines an "average" 
velocity within the pore. 

The continuity equation for one dimensional fluid flow in 
porous media with nonconstant porosity distribution is, 


D(po .) 
as ee 


Q9| a> 
vie 
Ul 
© 


, (tegen) 


Substituting in the Dupuit-Forcheimer assumption and Darcy's 


Law, the continuity equation becomes, 


E ile 9 (pe _) 
dP, tL a,lam_ il aydP_ iu a 
fee ho Cm ox) bx dx p.m 9t  ~ ° ‘TTt-22) 


Equation III.12 along with the boundary conditions can be 
integrated to provide the pressure and the pressure gradient 
distributions through the porous medium. The boundary condi- 


Bronms are P(0) = Poe the ambient pressure, and P(L) = PL: The 
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pressure differential across the porous medium is AP = PL pre | 
Knowing the pressure gradient distribution, the pore velocity 
distribution is obtained from Darcy's law and the Dupuit- 
Forcheimer assumption. Derivations of equation III.12 and 


the continuity equation are presented in Appendix A and B, 


respectively. 


C. SEMENOV MODEL OF COMBUSTION 

In this work we adopt the combustion model of N. N. Semenov 
as described in the texts of Frank-Kamenetskii [12] and Vulis 
[13]. A brief discussion of those features of the model which 
relate to the present investigation will be given. Fundamental 
to the model is the relation of reaction rate to temperature, 
and the interaction of heat generation and heat transfer. 
The reaction rate, Ras is represented by the Arrhenius law for 
a simple n-th order reaction, 


R = Ad" exp(—*) (III.13) 


¢ RT 

u 

where ATS is the characteristic time of the chemical reaction, 
E is the activation energy, RY is the universal gas constant, 
T is the absolute temperature, and > is the oxygen concentra- 
tion. A simple reaction is one in which the rate depends on 
the concentrations of the reactants and not on the products. 
The heat generated by the exothermal reaction is obtained by 


multiplying Ro by the heat of combustion. 
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In Figure III.2, heat generation Be is plotted versus 
temperature. The resulting curve is referred to as the S- 
curve due to its sinuous shape. The S-curves are distinguished 
by two regions. In region I, the reaction rate and tempera- 
ture are low and as a result, there iS an excess supply of 
oxygen. In this region the reaction 1s controlled by the 
temperature, and is called the kinetic regime of combustion. 
In the kinetic regime, the reaction rate increases exponen- 
tially with increasing temperature. In region II, the higher 
reaction rates are only slightly affected by the higher tem- 
peratures, and the reaction is limited by the availability of 
oxygen. Region II is called the diffusion regime of combus- 
tion. It should be noted that S-curves are not obtained by 
plotting RG versus temperature for a constant oxygen concen- 
tration. The flattening of the S-curve at the higher tempera- 
tures is due to decreasing oxygen concentration at those 
temperatures. The S-curve, which represents the behavior of 
an actual system, demonstrates the strong coupling between 
temperature and oxygen concentration. 

In addition to the S-curves, Figure III.2 shows two con- 
vection heat transfer curves, Fel and Too, at air flow tenm- 


peratures, T, and T respectively. The equation for the 


i Zn 


heat transfer lines is of the form, 


Geen tt. > 7.) (III.14) 
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where h is a convection heat transfer coefficient, and T. and 
T, are the solid surface and air temperatures, respectively. 
The intersection of the heat generation and heat transfer 
curves, RS = dos defines the stable quasi-stationary point A. 
For a system in thermal equilibrium, the solid and air tempera- 


tures will remain at TA and Th respectively. Perturbing the 


system from one equilibrium position to another causes the 
quasi-stationary point to adjust accordingly (e.g., raising 


the air flow temperature from T, to T will cause the quasi- 


Z 


Stationary point to move from T, to T At point I, the R 


A I): 


and GQ, curves are tangent, and an infinitesimal increase in 


J 


temperature will result in a jump in reaction temperature from 
T; to T- At point I, which is defined as the "critical igni- 
tion condition", the reaction moves from the kinetic regime 


Bemene diffusion regime of combustion. Temperature Tr is 
referred to as the ignition temperature. Although the des- 
cription of the process just presented is for a surface in 
which only heat transfer by convection occurs, the underlying 
ideas carry over for the total heat transfer at a point in the 
porous medium. The heat transfer can include contributions 
from conduction and radiation, as well as from convection. 

In his original work, Semenov [14] developed his theory 
to explain reaction behavior within homoegenous gas mixtures. 
Frank-Kamenetskii [12] later adapted this theory to describe 


the behavior associated with heterogeneous mixtures (i.e., a 


solid surface and air). In their study of coal combustion, 
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Thomas, Stevenson, and Evans [15] state that a "temperature 
jump", as described by Frank-Kamenetskil, may or may not 

occur depending on the partial pressure of oxygen. This pres- 
ent investigation does not consider this aspect of the com- 
bustion theory. However, analyses were made using the 

present model to determine whether the reaction takes place 

in the kinetic regime or in the diffusion regime. Moreover, 
the Semenov model will be used to explain these results in 


Chapter V. 


D. ARRHENIUS LAW OF REACTION RATE 
For carbon reacting in air, Parker and Hottel [16] pro- 
posed the following Arrhenius expression for the reaction 


rate in units of kg-carbon/m*-s, 


a) 


R = 9 SS Tae ut: exp (————) CT eS) 
C Wee ee wer : 
Gs u Os 





tj > 


where "Oy, 1s the partial pressure of oxygen in atmosphere 
units, and RA is the universal gas constant (1.986 cal/gmole- 
K). Equation III.15 assumes a simple first order reaction 

for the combustion of carbon in air. Frank-Kamenetskii [12] 
has shown that the experimental data of Parker and Hottel is 
better correlated by a fractional order reaction, where the 
reaction order, n, 1s between 1/3 and 2/3. Replotting Parker's 


and Hottel's data for n = 1/2, yields the reaction rate ex- 


pression used in the present analysis, 
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R = 2.065 x 10°(R, ayes eo) (iim. 16) 


- 2 R 
u ( 


In expression III.16, Ro romrnmUnits Of lbm-carbon/ft*- 
Ane "0. is the gas constant for oxygen (48.29 ft-lbf/lbm-R), 
® is in Lbm/ft>, RY is 1.986 Btu/lbmole-R, and P is in Rankines. 
The ideal gas law was used to replace the partial pressure, 
Bo. ° Figure III.3 and Table III.1 show the results of plotting 
Parker's and Hottel's data for values of n equal to l, 2/3, 
1/2, and 1/3. The values presented in Table III.1 are in S.I. 
bats. 

In order to determine the rate of heat generation and the 
Fate of oxygen consumption, the chemical reaction for the 
combustion process must be considered. For nonvolatile com- 


bustion of carbon and oxygen, two reactions that describe the 


process are, 
O > CO (Tite 17 > 
C+0 + CO (LET 2s) 


The ratio of the mass rates of carbon monoxide to carbon 
dioxide produced increases with increasing temperature. 
Arthur [17] presents an expression for the rate ratio as a 
function of temperature (in Kelvins). 


ASO Cogn es (III .19) 
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FIGURE III.3 --- Results of plotting reaction rate 


data of Parker and Hottel [16] for 
geveral values of n. 


82 


















mie a 
eee]. eS 
m-atm-se mole 


119 x 10° 27500 


sof 





TABLE III.1 --- Results of plotting reaction rate 
data of Parker and Hottel [16] for 


several values of n. 
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As a result of this dependence on temperature, the stoichio- 
metric ratio and the heat of reaction will also be functions 
of temperature. Defining the fraction of carbon monoxide 


being produced by, 


(CO/CO,) 
Foo = TF 1CO7COp (III.20) 


and the fraction of carbon dioxide as, 


sl 


160), T + (CO/CO5) CLiL 2 ds) 
the heat of combustion is then expressed as, 
AH, = Pao AHao + “CO, paxcoe CEEL «22) 
Values for the heats of combustion, AHao and AHao , as func- 


2 
tions of temperature were obtained from the JANAF (Joint Army, 


Navy, and Air Force) Tables [18]. The stoichiometric ratio 


(fuel to oxygen) of the overall reaction is, 


fR = uo “co, ‘feo, Fao + £a0 Fao ) (ia i.2 3) 


where Luo 1s the stoichiometric ratio for reaction III.17 and 


EO is the stoichiometric ratio for reaction III.18. The 
2 
rate of heat generation, nee and the rate of oxygen consumption 


can now be expressed by, 
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R = AH, R rni. 24) 
R = £ R EE 2 5) 


It should be pointed out that Parker's and Hottel's work [16] 
and Arthur's work [17] were accomplished using specific types 
of carbon. Smoot and Pratt [19] and Frank-Kamenetskii [12] 
list tables and references for many other rate expressions, 
depending on the particular type of carbon used. The present 
model was developed for any rate expression of Arrhenius type 
for carbon consumption. For this investigation, we will re- 
main with Parker's and Hottel's expression III.16 as modified 
by Frank-Kamenetskii. 

As previously stated, the particle diameter decreases as 
combustion progresses. The rate of decrease depends on the 
amount of carbon consumed at a point over time. Past observa- 
tions have shown that the effect is significant when the reac- 
tion 1s concentrated in a small region of the porous medium. 
To take this into account, expressions for the time rate of 
change of the diameter as a function of reaction rate were 


derived. The equation for cylindrical fibers is, 
aaa 2 
qd = <-2 R,z2D FiAGaIp od) (III.26) 


and for spherical particles, 


Ae -2R,2D°/(n o, a*) (III. 27) 
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where Po is the bulk mass density of carbon. The diameter 


and the reaction rate are functions of both time and position. 


E. HEAT TRANSFER EQUATIONS FOR POROUS MEDIA 

In a previous investigation of porous media, Green and 
Perry [20] developed heat transfer equations for the solid 
and fluid phases. The model included the three basic mechan- 
isms, (1) physical movement of the fluid with its own heat 
content, (2) conduction of heat through the solid and fluid 
phases, and (3) convective heat transfer between the solid 
and the fluid phases. Radiation heat transfer was assumed 
to be negligible, and combustion was not a consideration. 
The differential equations for the solid and the fluid phases, 


obtained from energy considerations, 








gT me: 
pC. (l-P) aE = k(1-p) ~ + h2(T_ - T) (iii. 23) 
2 
aT. oT 3 T 
eye 5e (CUP PLC Ua S + k P et - h2a(T - TS) (hPl29) 


and were solved numerically by Green and Perry using the 

finite difference method. In equations III.28 and III.29, 

C 1s specific heat, h is the internal convection heat trans- 
fer coefficient, k is thermal conductivity, and z is the wetted 
Surface area per unit volume. The s and w subscripts refer 

to solid and fluid properties, respectively. Their model does 
not account for the temperature dependency of the properties. 


Other investigators such as Sundaresan [21], Riaz [22], 
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Mendelsohn [23], and Schneider [24] have presented a single 
equation approach to describe heat transfer in porous media, 
dT enh 32D 

[(1-p) cp tpe, po Jae = -Up c Paz + ie hea te 2) (ole 50) 
Equation III.30 can be derived from equations III.28 and 
III.29 by assuming the solid and fluid phase temperatures are 
equal. Green and Perry proposed that the equal temperature 
assumption was valid for systems satisfying the following 


condition, 


hz 1/2 “w 


A = ES : 


> .342 Cine 3h) 
where a is the thermal diffusivity of the fluid. 

In the present investigation, the heat transfer equations 
are generalized to also include (1) radiation, (2) internal 
combustion, (3) temperature dependency of properties, and 
(4) compressibility effects of the air. 

An energy balance on the carbon gives the heat transfer 
equation, 


auc aT 


9 = Se ee = = z: ec 
axl (1 p) ue ax ie T,) ¥ AGT i Be sec at 


Sm eae), 
where kK. 1s a pseudo thermal conductivity to account for 
radiation between particles. The derivation of equation 

III.32 is presented in Appendix A. The effective conductivity, 


Kae of the porous solid was proposed by Russel [25], 
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ee. 33) 


273 


(l-p' wjel 


where Ke and Ke are the bulk thermal conductivities of carbon 
and air, respectively, and p' = (l-p). Russel's expression, 
based on an electrical resistance analogy, can be used for 
the full range of porosity, from 0 to l. 

The pseudo thermal conductivity due to radiation, kee is 
obtained as follows. Treating the idealized geometry of the 
porous medium as a series of closely spaced walls, the net 


radiation heat flux between two of the walls is, 
gq = se (Tt = ne) aS Cia} 
ie IL 2 ‘ 


where ¢€ is the emisSivity, and o is the Stefan-Boltzman con- 


stant. By expanding qd in a Taylor series about T an 


Wis 
analogy to Fourier's law of heat transfer by conduction yields 


the effective radiation conductivity as, 


k= 4066 T2/(2-e) (III. 35) 
Expression III.35 is then used for kK. iarequatieon Lil.32. The 
details of the above derivation are presented in Appendix B. 
The application of expression III.35 to porous media is also 
proposed by Rohsenow and Hartnett [26]. In amore rigorous 
development, Whitaker [27] presents an alternative approach 


for treating radiation heat transfer in porous media. 
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Following the experimental work of Yoshida, Ramaswami, and 
Hougen [28], the internal convection heat transfer coefficient 
is given by the empirical correlation, 


hea 0.91 Re’ °°? 


-2/3 
(pc. G(cu/Kk.) en ] (eel 36) 
where ~ 1s equal to .91 for cylindrical fibers, and equal to 
1 for spherical particles, G is a pseudo mass velocity given 
by pe iu, and C. is the specific heat of air at constant pres- 
sure. The fm subscript refers to the properties evaluated 


at film temperature. Re' is a pseudo Reynolds number defined 


by, 
Re' = G/(zuy) (III. 37) 


The air properties vary with temperature, and h varies with 

temperature. Finally, the reaction rate for the heat genera- 

tion term in equation III.32 is given by expression III.24. 
An energy balance on the air within the porous medium 


provides the second heat transfer equation as, 





i 2) - ae ae ye ee 
ax P* 35x PP aM ax Cc a aox et 5 (ame 
Cire 3S } 


The density of air, Oo: 1S approximated in the model by the 
ideal gas law. The term, 3(pP)/3x, is due to the compressi- 


bility of the air. The derivation of equation III.38 is 
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presented in Appendix A. All properties in equation III. 38 

vary with temperature. The properties of standard air were 

used in the analysis. During combustion, oxygen in the air 

is replaced by carbon monoxide and carbon dioxide. Calcula- 
tions showed that the gas mixture would have thermophysical 

properties slightly different from those of standard air. 

In a “worst case" situation, the average difference between 

the thermophysical properties of air and a mixture of .79 No 


and .21 CO, is 7 percent. This mixture assumes that the oxy- 


2 
gen is totally consumed. Viscosity had the largest difference 
of 13 percent, and specific heat the smallest at 1.5 percent. 
The properties of a .79 N, and a .21 CO mixture are approxi- 
mately those of air (different by less than 2 percent). At 
typically observed temperatures (greater than 1200 deg-F), the 
combustion gas will be mostly carbon monoxide, and the error 
introduced by assuming the properties of standard air is mini- 
mum. A small difference in the properties is acceptable since 
an additional molecule mass transfer equation for either carbon 
monoxide or carbon dioxide would be necessary to calculate 

the effective properties of the gas mixture. The methods used 
to obtain the thermophysical properties of the gas mixtures 
above were those of Mason and Saxena [29], and Reynolds and 


Perkins [30]. The polynomial expressions used to calculate 


the thermophysical properties of air are derived in Appendix B. 


FF. OXYGEN DIFFUSION EQUATION FOR POROUS MEDIA 
As a result of combustion, the heat transfer equations 


include four response variables, the carbon and air temperatures, 
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T. and Toe the concentration of oxygen, ¢, and the total 
internal pressure, P. The fourth field equation necessary 
to complete the system is obtained from species conservation 
considerations. The oxygen molecule transport mechanisms 
included in the model are (1) molecular diffusion resulting 
from concentration gradients (Fick's Law), (2) convective 
mass flow, and (3) oxygen consumption due to combustion. 
Diffusion resulting from pressure and temperature gradients 
was considered negligible. An example cited by Bird, Stewart, 
and Lightfoot [21] where pressure diffusion is important is 
the centrifuge separator; and for thermal diffusion, the 
Clusius-Dickel column, where very steep temperature gradients 
are used to separate complex mixtures of organic molecules. 
Thus, the oxygen molecule transfer equation for a porous 
medium is, 

0 .) 


dd) Oo 2 = oo 
S-(p 0, $2) ay (PU >) Ry 2 = P x¢ (Pers 9) 


Q 


The derivation of equation III.39 is presented in Appendix A. 
The second order term in equation III.39 results from Fick's 
law of diffusion. The effective diffusivity of the porous 
medium, Dos is a function of pressure, temperature, and pore 
geometry of the medium. Scheidegger [9] presents several 
models for ee which have been proposed by a number of inves- 
tigators. Briefly, there are two limiting cases of diffusion, 
(1) diffusion associated with the collision between gas mole- 


cules, and (2) diffusion associated with the collision of 
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gas molecules with the pore walls of the medium. This second 
phenomenon is discussed by Bennett and Myers [32], and is 
referred to as Knudsen diffusion. The former case occurs 
when the mean free path of the molecules is smaller than the 
pore diameter. The expression for mean free path, w, given 


by Treybal [33] is, 
WwW = ac ae (R Bay ira M) 1/2 (ITIIT.40) 


where J. 1s the gravitational constant and M is the molecular 
weight. From equation III.40, a value of w = 6/100 is obtained 
for oxygen for the geometry of typical porous media encoun- 
tered in the analysis. Thus, collision between molecules is 
the governing diffusion mechanism. A semi-empirical expres- 
Sion, proposed by Gilliland [34], is used to obtain the diffu- 
een coefficient of oxygen into air. Gilliland's expression 

1S given by, 


D9 = 435.7 7°73 (M_* +MO*) /{P (V2 /3 + v0) 1 (III. 41) 
a a 0. 


25 


where 90 1S in units of em“/sec, P is the total pressure in 


Pa, ‘a and Vo are the molecular volumes of air and oxygen, 
2 


respectively. M. and Mo are the molecular weights of air and 
2 


oxygen. The values of Mis and Vo were obtained from Holman 
2D 
[35] as 29.9 and 7.4, respectively. The effective diffusivity 


proposed by Denbigh and Turner [36] for a porous medium is, 
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D, = 0/1 (III. 42) 
Expression III.42 accounts for the tortuous path the oxygen 
molecules follow through the porous medium. Lastly, the oxy- 
gen consumption term in equation III.39 is given by expression 


re. 25. 


G. THE SURFACE RECESSION PROBLEM 

The previous system of equations describes the combustion 
process occurring within the porous medium. As combustion 
progresses, the oxygen concentration goes to zero, and as a 
result, the reaction moves to the air inlet surface of the 
porous medium. Thereafter, the carbon consumption takes place 
at the surface. During the surface recession phase, the thick- 
ness of the porous medium decreases with time, 1.e., surface 
recession. This suggests partitioning the problem into two 
parts, (1) initial combustion within the porous medium, and 
(2) combustion at the air inlet surface. The second phase 
of the problem is formulated as a moving boundary problem. 
The governing equations for the moving boundary problem are 
those developed previously with the exception of the oxygen 
molecule mass transfer equation. With combustion occurring 
at the surface, oxygen is totally consumed in a shallow region 
near the surface. Experimentation by Koizumi [37] has shown 
this penetration of the oxygen to be 1 to 2 particle diameters. 
As a result, the oxygen molecule mass transfer equation is 
not needed for this phase of the problem, and the heat genera- 


tion can be treated as a planar source. 
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The carbon heat transfer equation with the modified heat 
generation term is, 
oT 


-20 S oo! 
L aah (l-p) Wea eae -hz (T.-T.) inl eae (n = 0) 


oT 
oy 


notes 
= cme ecca Lan ot 


VEEL 43) 
where n is equal to x/L. The first term on the right side of 
equation III.43 arises from the x coordinate becoming a func- 
tion of time during surface recession. A similar term appears 
in the air temperature and combined Darcy's law-continuity 
equations as well. Transformation of the field equations from 
a fixed coordinate to a moving coordinate system is presented 
in Appendix A. As was observed (see Figure V.2), the tempera- 
ture gradient through the porous medium during surface reces- 
sion was small. As a result, the additional term may be 
omitted without affecting the results. The moving boundary 
problem formulation is based on Crank's [38] extension of a 
method proposed by Murray and Landis [39]. The thickness, L, 
is updated continuously by evaluating its time rate of change 


given by, 


L = -(puf,) pe (a) OQ] (III.44) 
x=0 x=0 


Expression III.44 is obtained from amass balance on the 


Carbon (1.e., the time rate of change in carbon is equal to 
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the carbon consumption). Noting that all the oxygen entering 
the plate is consumed in a small region near x = 0, the carbon 
consumption can be represented in terms of the oxygen flow 

and the stoichiometric ratio (l.e., Re = -puf,9,A where A 

is area). The time rate of change for carbon is (1-p) p AL. 
The planar heat source, Hoo in equation III.43, is determined 
by the amount of oxygen entering the porous medium at x = 0. 
This is converted to heat generation by using the stoichio- 
Thus in 


metric ratio, £ and the heat of combustion, AH 


Re R° 
equation III.42, the planar heat source becomes put ,AH,¢.- 
Ozisik [40] presents a general discussion and references for 
the moving heat source problem. 

Transition from the internal combustion problem to the 
Surface recession problem occurs when the porosity at x = 0 
nears a value of 1. The computer program, implementing the 
analysis, was written so that the change from an internal 


combustion problem to a surface recession problem occurs 


automatically. 


H. BOUNDARY CONDITIONS 

Three sets of boundary conditions may be used for equations 
Mee s2, LIL.38, and 11.39. Each set of boundary conditions 
approximates a physical situation (e.g., thermal flow reactor 
Or Fontenot's experiments [1]). The boundary conditions 
were the nearest approximations that could be made and remain 
within the limitations of a one dimensional model. MThe first 


and second set of boundary conditions are typical of thermal 
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flow reactors with and without radiation from 


surfaces. These are 


and 


aT 


- Cc. u(T. bes) 


Q 
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with radiation, 
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the boundary 


= 0 (TEpe4 5) 
= L (III.46) 
= 0 (lita) 
= L Ciao) 
= 0 (TET .49) 
= L (III.50) 
= 0 (Tole. ) 
= L (rT . 52) 
= 0 ey Ene) e 
= L (ET. 54 } 
= 0 elie 5) 





r= 0 et (III. 56) 


Expressions III.45 and III.46 represent insulated boundaries 

on the porous solid (1i.e., no heat loss from the porous solid 
to the environment). Expressions III.5l1 and III.52 provide 

for radiation heat transfer from the porous solid boundaries 

to the environment. The boundary conditions for the air 
temperature and oxygen concentration, expressions III.47 
through III.49 and III.52 through III.55, are the Danckwerts' 
boundary conditions for flow reactors. A convincing discussion 
of the Danckwerts' boundary conditions applied to mass diffu- 
Sion equations for porous media is given by Bischoff [41]. 

A brief summary of the discussion was presented by Vatikiotis 
[2]. An analogy to Bischoff's approach for the air heat trans- 
fer equation is presented in Appendix B. 

The third set of boundary conditions approximates the situa- 
tion for Fontenot's experiments [l]. In the experiments, plate 
specimens were mounted into the wall of a wind tunnel. The 
interior surface of the plate (i.e., at x = L) was exposed to 
the wind tunnel air flow, and the external plate surface (l.e., 
at x = 0) was exposed to the outside environment. The expres- 
Sion for calculating the pressure differential across the 
porous medium as a result of the external flow over the surface 
at x = L is formulated in Appendix B. It must be emphasized 
that the purpose for modelling Fontenot's experiments was to 
analyze the combustion behavior within the plate specimens 


using the one dimensional model. However, the external flow 
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and plate surfaces involved convection heat transfer requiring 
two dimensional analysis. Noting this, the following approach 
was adopted to provide for convection heat transfer at the 
boundaries of the one dimensional model in a qualitative 


sense. The boundary conditions are, 


qT - 
oa - = 4 ” 

(kK tk Jae = hy (7-7) to¢ (P_-T,) x = 0 (III.57) 

ae A4 r4 
(kK +k) = -h, (T,-T,)-9¢ (T-T,) x= L (III.58) 
a = T,, <= 70 (III.59) 
oT. 
= (OO x = L (III.60) 
—- >. x = 0 (elr 61) 
-: Be 56 x=L (III .62) 


As stated above, the values of hy and h. in the convection 
terms of expressions III.57 and III.58 are difficult to inter- 
pret for a one dimensional model. The difficulty arises from 
considering the boundary layers on the external surfaces of 
the porous medium (e.g., considering the porous medium as a 
flat porous plate, the boundary layer resulting from natural 
convection on the x = 0 surface increases in thickness in the 
vertical direction, and for forced convection at the x = L 


surface, the boundary layer thickness increases in the flow 
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direction). For this investigation, Ay was approximated 
by, 
Ka ligO 
By eae (een .63) 
x=0 

This expression is based on the analytical results presented 
by Merkin [42]. Expression III.63 is obtained by considering 
natural convection heat transfer from a vertical flat plate 
with "suction". From Merkin's results, the boundary layer 
thickness becomes constant at some distance in the vertical 
direction, and thus, the convection heat transfer may also be 
considered constant. The magnitude of the convection heat 
transfer coefficient, Ao, varies in a direction parallel to 
the external flow, UL. at the x = L surface. In addition, hao 
depends on the efflux of the gas at the surface. To simplify 
the analysis, he was approximated by the relations for a smooth 


flat plate given by Holman [35] as, 


LVS: 1,7 


k 
h, = .664 —2 pr Re (III.64) 
for laminar flow, where L is the distance from the boundary 


Payer initiation, and, 


a ibys 8 


(3037 Re = 850) Cir 265) 


for turbulent flow, where Pr is the Prandtl number. The 


Reynolds number here is based on the external flow velocity 
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over the surface at x = L. Kays [43] provides an alternate 
scheme for treating boundary layers with "blowing and suction". 
The air temperature and oxygen concentration boundary condi- 
tions, III.59 through III.62 were also selected because of 

the difficulties in treating external surface boundary layers 
in a one dimensional model. Alternatives to the essential 
boundary conditions of the air and oxygen are the Danckwerts' 
conditions. However, using the Danckwerts' conditions with 
convection boundary conditions poses the problem of approxi- 
Mating a reasonable value for u at x = 0 and x = L (i.e., 
boundary layers affect the component of velocity, u, normal 

to an external surface). To overcome these difficulties, it 
seemed prudent not to have pore velocity and the convection 
heat transfer coefficients, Ay and ho, appear in the same set 
of boundary conditions. However, the difficulty still re- 
mains of specifying a distance from the origin of the boundary 
layer for expressions III.64 and III.65. In the initial work 
by Vatikiotis [2], an arbitrary distance of 1 foot was used 


to compare results. 


fer LNITIAL CONDITIONS 

To initiate combustion, the porous medium must be brought 
to a temperature at which the heat generation rate is suffi- 
ciently high. Experimentally, there are many techniques for 
doing this. However, it is difficult to measure the actual 
temperature and oxygen concentration during the experiments. 


This difficulty is carried over to identifying a reasonable 
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set of initial conditions. Noting this, the model was de- 
veloped so a problem can be started at ambient conditions 
with a heat flux applied to the carbon at a boundary. The 
heat flux is mathematically treated as, 


oT 


(1-p) (k, + k)) = = -G_. Elie) 


This approach of starting the problem simplifies the task of 
determining a reasonable set of initial conditions. In addi- 


tion, any arbitrary set of initial conditions can be specified. 
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IV. IMPLEMENTATION OF NUMERICAL METHODS 


A. SOLUTION ALGORITHM 
The field equations and auxiliary equations which define 
the problem were presented in the previous sections. The 


following is the scheme used to solve for the dependent 


variables Toe T 


i. 


~ Oo Wee, Ly, and p.: 

Starting with the Dupuit~-Forcheimer expression III.10, 
the filter velocity, Q, 1s represented as a function 
of the pore velocity, u. 

Substitute for Q from the previous step into Darcy's 
law, equation A.l. 

The expression formed in step 2 is solved for the pore 
velocity, u, as a function of dP/dx. This procedure 
yields equation A.2. 

Equation A.2 1s substituted into the continuity equation 
A.3 resulting in the 2nd order partial differential 
equation A.4 with the pressure, P, and the air density, 


fe as the dependent variables. 


a? 
Expanding equation A.4 yields equation A.5. Equation 
A.5 1s one of 4 equations which were integrated to 
produce the problem solution. By taking the air den- 
cy , Pp,, as a constant over a time step of integration, 
equation A.5 becomes an ordinary differential equation 


with respect to x. For equation A.5, the shooting method 


(with Euler's formula for the integration algorithm) 


52 





0. 


was used to solve for the pressure, P(x), and the 
pressure gradient, dP(x)/dx at each time step. The 


air density, p_, was updated at each time step using 


a 
the ideal gas law, equation B.4l. The rate of change 
of the air density was neglected during the first two 
time steps of integration. Reasons for neglecting the 
time rate of change of the air density are given in 
Appendix C.2. 

Having solved equation A.5 for the pressure gradient, 
dP/dx, the pore velocity can be obtained from equation 
Ac 2. 

In solving the air heat transfer equation III.38 and 
the oxygen molecule transfer equation III.39, the pore 
velocity, u, and the air density, P,, are obtained from 
equations A.2 and B.4l, respectively. Equation III.38 
and III.39, along with the carbon particle equation 
III.32, were transformed to ordinary differential equa- 
tions in time by a Galerkin FEM formulation. 

The ordinary differential equations from step 7 are 
given by equations C.17 through C.19. The time inte- 
gration was performed by a version of Gear's method 
developed by Franke [44]. 


All thermophysical properties (i.e., k k oe 


k 
eo fs a a 


c u) being functions of temperature are continuously 


ca’ 
updated during the transient analysis. 
As carbon is consumed during combustion, the porous 


medium properties (1.e., p, m, 46, dad) change with time 


and are updated continuously as well. 
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12. No attempt was made to match the orders of convergence 
of the integration algorithms. 

Implementing the above solution algorithm resulted in a 
computer program for the solution of the problem. Several 
options were incorporated in the computer program allowing 
the user flexibility in the analysis. Boundary conditions 
Simulating fixed air inlet temperature, and convection and 
radiation heat transfer from the exterior surfaces can be se- 
lected. The problem can be initiated in two ways; an arbitrary 
set of initial conditions can be specified, or a heat flux, 
shown by expression III.66, can be specified with the porous 
medium at ambient conditions. In addition, the porosity at 
x = 0 causing transition to the surface recession phase can 
be specified. Other options permit investigating materials 
other than carbon (e.g., boron), studying the effects of non- 
constant distributions of porous media properties. 

The computer program includes a main program and 24 sub- 
routines (7 of the subroutines comprise Franke's [44] inte- 
gration routine). The function and description of each of 
the subroutines are found in the computer program listing. 
Figure IV.1 shows the flow of the program, and briefly des- 
cribes the process at each step. Each block represents a 
subroutine. The configuration of the blocks (1i.e., blocks 
within blocks) in Figure IV.1 show the six levels of the 
computer program. In addition, other processes associated 


with the integration routine occur inside the block shown by 
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Main program. 








Read in parameters, initialize flags and 
counters; Calculate initial properties 
and pore velocit 


Construct FEM nodal point/element 

correspondence array. 

Print current response variables 

and properties. 

Check for transition TRUE Reset flags. counters delete 
to surface recession. @ (called once). 


FALSE 










Call integration routine for T., T, P : CD deleted 
c* a 
during surface recession). 


Calculate properties. 


Test for oxygen instabiltiy. 


If surface recession, update thickness, L. 
Calculate pore velocity and pressure gradient. 


If surface recession, transform reaction rate 


expression to surface flux. 





Set flag for surface recession if p(n=1l) = pmax. 


Zero-out system matrices. 
Generate FEM operators. 


Form mass and system matrices. 


Adjust matrices for boundary conditions. 





. * . 


Adjust matrices for oxygen instability. 


| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 
| 






Perform integration. 


FALSE— Test flags for program termination. 


TRUE 






Print last values of response variables 
and properties. 





FIGURE IV.1 - Flow 
of computer program. 
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the dashed line. These processes are discussed in detail in 


Franke's report. 


B. STIFFNESS CONSIDERATIONS 

Initially, the ordinary differential equations resulting 
from the finite element formulation were integrated explicitly 
using a sixth order Runge-Kutta method (IMSL subroutine DVERK). 
The largest obtainable time step was approximately 1a 
seconds. A larger time step caused the integration to become 
unstable. Gear [45] described this behavior as typical of 
"stiff" differential equations, and suggested that implicit 
integration methods be used to improve the efficiency of the 
integration. Gear's method, presented by Brown and Gear [46] 
and modified by Franke [44], was adopted as the integration 
method. As a result, time steps of several seconds were ob- 
tained for some problems (i.e., for surface recession and when 
extinction occurred). In addition, moving the exponential 
reaction rate from the excitation vector to the stiffness 
matrix (discussed in Appendix C.1) improved the integration 
efficiency. The improved efficiency was attributed to a better 
approximation of the Jacobian matrix in Gear's method. 

Bui [47] and Burka [48] discuss the difficulties of solving 
differential equations exhibiting "stiffness" and propose 
alternatives to Gear's method. For a system exhibiting "stiff- 
ness", small time steps are needed for the rapidly responding 
terms while the integration must be continued for a long period 


to account for the slowly responding terms. A measure of 
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"stiffness" for a system of differential equations is the 
ratio of the real part of the maximum eigenvalue of the Jacob- 
ian matrix to the real part of the minimum eigenvalue. This 


is called the “local stiffness ratio", and is represented 


by, 
Max JA. | 
. Peiiynes 7 n 
S = Sawai = fA. | (ce 1 ) 
Metl pcee)5 11 - 


where n is the number of equations. A system of equations can 
be considered "stiff" if S has a value greater than 10, and 
the real parts of the eigenvalues are negative. 

Lastly, most numerical integration methods are bound by 
a maximum time step for which the solution remains stable. 


This is defined by, 


JAt-A,| < kK (IV.2) 


where At is the time step, As is the real part of the maximum 
eigenvalue of the Jacogian matrix, and K is a constant depend- 
ing on the integration method (usually 1 < K < 10). For 
example, Bui [47] states that Euler's method requires 

JAtA, | eee oe re rore, for a maximum Jacobian matrix eigen- 
value (representing fastest response) of 50 (seconds) +, the 
largest timestep for the integration is .02 seconds. Exceed- 
ing this value results in the integration becoming unstable. 


This explains the maximum timestep of Vor seconds for the 


Sixth order Runge-Kutta method as previously discussed. 


5) 7 





C. TREATMENT OF NUMERICAL DIFFICULTIES ASSOCIATED WITH 
OXYGEN CONCENTRATION 


In problems of combustion or kinetics, where dependent 
variables go to zero, numerical difficulties arise. In his 
discussion, Frank-Kamenetskii [12] points out that for frac- 
tional order reactions (expression III.16), both the oxygen 
concentration and concentration gradient can become equal to 
zero within the porous medium. This requirement is difficult 
to satisfy with approximate solution methods. In the initial 
effort (Vatikiotis [2]), run time errors resulted when the 
oxygen concentration became slightly negative. This was 
caused by attempting to take the logarithm of a negative con- 
centration in the fractional order term in expression III.16. 
One approach for correcting this problem was to use the abso- 
lute value of the oxygen concentration in the reaction rate 
term. This resulted in the concentration becoming increasingly 
negative. Another approach was to check the values of oxygen 
concentration at each integration interval and setting the 
concentration to zero for the negative values. At the next 
integration, the reaction rate was zero at the locations of 
zero concentration. Without consumption the oxygen concentra- 
tion at the nodal point previously set to zero becomes positive 
and the oscillation repeats itself. It was observed that de- 
creasing the integration time step and the length of the ele- 
ments (i.e., increasing the number of nodal points) minimized 
the oscillations. Since CPU time increased for a given problen, 


the approach was not totally satisfactory. However, this 
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method was adopted for obtaining the results presented by 
Waeikiotis (2). 

An attempt to resolve the above numerical difficulty for 
the present model was by implementing the Moving Finite Element 
method of Gelinas, Doss, and Miller [49]. In the method, 
the nodal points defining the elements converge to regions 
where they are needed to minimize the error in the numerical 
solution. As a result, the size of the elements in a region 
of large reaction rate become small, and move with this region 
in time. The motivation for implementing the Moving Finite 
Element method was Frank-Kamenetskii's [12] discussion of 
fractional order reactions as stated above. Also, it was ob- 
served that the oxygen concentration approached zero ina 
small interval within the porous medium. This behavior resulted 
in steep concentration gradients. It was thought that a num- 
ber of small elements in this region would improve the accuracy 
of the solution. The results obtained from the Moving Finite 
Elements showed the method was not effective in resolving the 
numerical difficulties discussed above. The method was subse- 
quently abandoned because of the additional expense associated 
with solving the differential equation governing the nodal 
Beant locations. 

The method adopted in the present model for minimizing 
the difficulties with the oxygen concentration is as follows. 
The time derivative of the oxygen concentration at a nodal 
point is set to zero when the concentration reaches a speci- 


6 


fied positive minimum value (less than 10. lbm/£t?) . At 
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the next integration interval, the time derivative remains 
zero if the concentration at the adjacent nodal point in the 
negative direction (1.e., upstream) has decreased. If the 
concentration at the adjacent nodal point has increased, then 
the time derivative is released. This method, in conjunction 
with Franke's [44] integration routine, appears to act in 

an iterative manner continually improving the solution. The 
values of concentration in the regions where the oxygen has 
been totally consumed are typically on the order of ome 
lbm/ft>. Although they do not give the details, Scheisser 


and Stein [50] have used the above method in their work of 


Simulating coal conversion. 


D. AN ALTERNATE SOLUTION STRATEGY 

The experience and difficulties of implementing the numeri- 
cal methods, as discussed above, have provided insight for an 
alternate strategy for solving the combustion and heat trans- 
fer problem. The alternate strategy consists of formulating 
the problem with two moving boundaries. One moving boundary, 
as in the present formulation, would account for surface 
recession (i.e., the porous medium thickness changing with 
time). The second moving boundary would account for the 
oxygen concentration going to zero within the porous medium. 
This second boundary would be located where the concentration 
becomes zero, between x = 0 and the first moving boundary. 
The advantages of the alternate solution strategy are (1) 


eliminating the numerical difficulties associated with the 
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oxygen concentration (discussed in the previous section), 

(2) relaxing the planar heat source assumption for the surface 
recession phase (discussed in Section III.G), and (3) possibly 
improving the numerical "stiffness" characteristics of the 


integration (discussed in Section IV.B). 
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V. RESULTS AND OBSERVATIONS 


fee BASIS FOR THE RESULTS 

The results presented in the following sections are based 
on the idealized porous medium described in Section III.A and 
Figure III.1l. For the idealized porous medium, porosity and 
permeability were explicit functions of pore diameter and unit 
cell thickness. Other porous media will have different geo- 
metric relationships for describing the properties. In addi- 
tion, for naturally occurring porous media, experimental methods 
would be employed to determine the geometric properties. 
Noting this, the results of this investigation are intended 
to describe, ina qualitative sense, the behavior of a carbon 
porous medium during combustion. As will be shown, the tempera- 
tures observed during combustion are highly dependent on the 
geometric properties and the pore velocity. Although, the 
temperatures are typical of those reported for carbon combustion 
experiments, it would not be reasonable to perform a quanti- 
tative comparison between the results obtained by the model 
and those obtained experimentally unless the porous medium 
properties were the same for both. This point will be demon- 
Strated in the next section. 
B. COMPARISON OF COMBUSTION MODEL RESULTS TO EXPERIMENTAL 

RESULTS 

The behavior of the temperature and oxygen concentration 


of the combustion and heat transfer model is in fair agreement 
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with the behavior reported in the experimental investigations 
@eeVulis [13], Koizumi [37], Fontenot [51], Spalding [(52}, and 
Kolodstev [53]. A quantitative comparison of the results can 
not be made for two reasons. First, there have been few 
experimental investigations which present data for transient 
behavior of combustion in porous media. Most of the analyses 
present quasi-steady data (quasi in the sense that for some 
period of time, there is little change in the response varia- 
bles). When transient behavior is reported, it is usually 
presented ina narrative form. Secondly, the reported des- 
criptions of the porous media used in the experiments are not 
sufficient. The properties of the carbon used (e.g., thermal 
conductivity, density, Arrhenius rate law, etc.) and the 
porosity and permeability relationships are needed to accurately 
Simulate the experiment. The profiles of the quasi-steady 
temperatures obtained by the combustion model and those ob- 
tained by Kolodstev's experiments will be shown together. The 
purpose of this is to demonstrate that although Kolodstev's 
porous medium could only be roughly approximated, the tempera- 
ture profile obtained by the combustion model was similar in 
shape and magnitude to that obtained experimentally. 

In the experimental investigations of Vulis [13], Fontenot 
[51], and Spalding [52], it was observed that for a given 
temperature, sustained combustion occurred for a specific range 
of velocities. Velocities outside the range, either greater 


Or lower, would cause the combustion to extinguish and the 
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Borous medium to cool to ambient temperature. This phenomenon 
was also observed for the combustion model. The explanation 
of this behavior, discussed in Section V.G, 1S summarized as 
follows. As pore velocity increases, convection heat transfer 
and oxygen supply for combustion both increase. Similarly, 

a decrease in pore velocity will cause a decrease in convec- 
tion heat transfer and oxygen supply. For a given set of 
conditions, pore velocities above the range will result in 

the heat transfer or heat loss overcoming the heat generation. 
Also, for pore velocities below the range, the oxygen supply 
for combustion is decreased to the extent that the heat genera- 
tion is overcome by the heat loss. In other words, the higher 
pore velocities tend to "blowout" the reaction (e.g., blowing 
out a candle), and the lower pore velocities tend to "starve" 
the reaction (e.g., placing a container over a candle). 
Important at lower pore velocities are the effects of heat 
transfer from the external surface of the porous solid (l.e., 
at the boundaries). The experiments performed by Vulis [13], 
Fontenot [51], and Spalding [52] allowed heat transfer from 
the surfaces by either convection or radiation. The effects 
of heat transfer at the boundaries on combustion are discussed 
mo oection V.G. 

Kolodstev [53] investigated combustion gas dynamics ina 
porous medium comprised of spherically shaped carbon particles. 
This is the same type of carbon used by Parker and Hottel [16]. 
In his results, Kolodstev presents the quasi-steady air tempera- 


ture as a function of the depth of the porous medium (particle 
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temperature was not reported). Figure V.1 shows Kolodstev's 
results plotted with the results obtained by the combustion 
model. Kolodstev's experiment could not be simulated accu- 
rately. As stated previously, the description of Kolodstev's 
porous medium was not sufficient. The description of the 
porous medium and the ambient conditions used to simulate the 
experiment is shown in Table V.1. Also shown are the known 
parameters given by Kolodstev. In Figure V.1, the temperature 
for the first inch of the 7.5 inches of the porous medium are 
shown (temperatures at greater depths were not reported). As 
can be seen, there is fair agreement in the shapes of the 
profiles. In addition, the magnitude of the temperatures pro- 
duced by the combustion model are similar to those obtained 

in the experiment. 

The combustion model results showed that once the reaction 
moved to the air inlet surface of the porous medium, the oxy- 
gen concentration penetrated the porous medium by 1 or 2 
particle diameters. This observation was also reported in 
the experimental results of Koizumi [37] and Kolodstev [53]. 
The response of the oxygen concentration, in general, will be 


discussed in the next section. 


Ge EXAMPLE RESULTS 

The intent here is simply to demonstrate the operation of 
the computer program. The behavior of a single porous medium 
subjected to several initial conditions is presented. The 


description of the porous medium and the ambient conditions 
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TABLE V.1 --- Geometry of porous medium and ambient conditions 
for simulating Kolodstev’s (€53] experiment. 


eParticle Shape: spherical 
@Particle Diameter (d): ~126 in 

Unit Cell Thickness (0): .126 in 
eSpatial Thickness of Porous Medium (LL): Pao in 

Porosity (Pp) 476 

Permeability (m) a Oe 19° 5° 

Bulk Thermal Conductivity of Carbon (k,): 86.0 Btu/ft-hr-F 

Bulk Specific Heat of Carbon (C,): ecxol Btu/lbm-F 

Bulk Density of Carbon (DO): 70.5 lbm/#t> 

Thermal Emissivity of Particles (€): ae 

Ambient Temperature (7],): 80.0 deg-F 
@Ambient Pressure (P)): 14.7 psi 
@Pressure Differential (AP): -.29 psi 
@Ambient Oxygen Concentration (Dy): OL Lbm/tt> 


® Known parameter for Kolodstev’s experiment. 
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used for the following problems are shown in Table V.2. The 
boundary conditions for the examples in this section were 
insulated boundaries on the porous solid temperature (ex- 
pressions III.45 and III.46), and the Dankwerts conditions on 
the air temperature and oxygen concentration (expressions 
meme 4? through III.50). 

As the first example (referred to as example A), a heat 
mux of 3.0x 104 Btu/ft?-hr was applied at x = 0 for 15 seconds. 
The porous medium was initially at a constant temperature of 
80 degrees Fahrenheit and a someten oxygen concentration of 


Seig2 lbm/ft> 


(ambient conditions). The results of this prob- 
lem are shown in Figures V.2-V.4. Figure V.2 shows the tem- 
perature increasing after the heat flux was removed. At 25 
seconds, the porosity at x = 0 reached .95 and the problem 
transitioned to the surface recession phase. The values of 
thickness for the times shown in Figure V.2 are .25 inches to 
25 seconds, .226 inches at 72 seconds, .133 inches at 357 
seconds, and .0238 inches at 517 seconds. The computer pro- 
gram was written such that the problem terminates when the 
thickness becomes 10 percent or less of the initial value. 
The air temperature shown by the dashed line in Figure V.2 
follows the carbon temperature with the exception of the air 
inlet region at x = 0. Figure V.3 shows the oxygen behavior 
for this problem. At 25 seconds, the oxygen concentration at 


: lbm/£t> and the penetration of the oxygen 


x = 0 was 3.6x10— 
was to x/L = .0l. This supports the assumption of using a 


planar source for the heat generation and deleting the oxygen 
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TABLE V.2 --- Geometry of porous medium and ambient conditions 
for the example problems in Section V.C. 


Particle Shape: 

Particle Diameter (d): 

Unit Cell Thickness (0D): 

Spatial Thickness of Porous Medium (L): 
Porosity (p>) 

Permeability (™M™) 

Bulk Thermal Conductivity of Carbon (k.): 
Bulk Specific Heat of Carbon (C,): 

Bulk Density of Carbon S/OR 

Thermal Emissivity of Particles (€): 
Ambient Temperature (7 ): 

Ambient Pressure (P)): 

Pressure Differential (AP): 


Ambient Oxygen Concentration (D,) 3 
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concentration equation during the surface recession phase. 
Figure V.4 shows the decrease in the thickness of the porous 
medium with time. The nondimensional parameter, A, given by 
expression III.3l, had an average value of .63 at the start of 
the problem. This value of A 1s above the minimum value of 
.342 proposed by Green and Perry [20] when assuming equal 
carbon and air temperatures. This suggests that a single heat 
transfer equation would have been sufficient to describe the 
temperature response for this problem. However, this is not 
valid for all cases. For example, when d, D, and L are 
doubled, that is, .0l inches, .0l1 inches, and .5 inches, 
respectively, A equals .22, and for d, D, and L equal to .025 
inches, .025 inches and 1.25 inches, respectively, \X equals 
-056. Hence, for the general case, the carbon and the air 
temperatures should be treated as separate response variables. 
As a second example (referred to as example B), a heat 
pax Of a0 % ave Btu/ft~-hr was applied at x = 0 for 14 seconds. 
The porous medium was initially at the same ambient tempera- 
ture and oxygen concentration as in the previous example. 
Figures V.5 and V.6 show the results of this problem. The 
porous medium cooled to ambient conditions in approximately 
41 seconds after the heat flux was removed. Applying the 
heat flux for 14 seconds was not sufficient to produce the 
conditions needed to sustain combustion. However, increasing 
the duration of the heat flux for an additional second results 
in sustained combustion. Figures V.5 and V.6 show the results 


starting at 12 seconds. The behavior of the system prior to 
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12 seconds is the same as in Figures V.2 and V.3 of the pre- 
vious problem. 

As the last example (referred to as example C), the problem 
was started with the porous medium at constant temperature 
and oxygen concentration of 1100 degrees Fahrenheit and .006 
lbm/f£t?, respectively. Figure V.7 shows the large cooling 
effect of the air entering the porous medium at x = 0 during 
the early part of the problem. As a result of this cooling, 
the heat generation began to raise the temperature interior 
to the porous medium before moving as a wave to the air inlet 
surface. The problem was terminated at 49 seconds with a 
porosity at x = 0 of .96. This was at the point of transition 
to the surface recession phase. Figure V.8 shows the behavior 
of the oxygen for this problem. At 49 seconds, the oxygen 


4 ipm/ft> and the penetra- 


concentration at x = 0 was 4.1x10 
tion of the oxygen was to x/L = .01. This again supports the 
assumption of treating the heat generation as a planar source 
and deleting the oxygen concentration equation during the 
surface recession phase. 

Other data generated by the model are spatial distributions 
of the following internal properties: particle diameter (d), 
porosity (p), specific permeability (m), pressure (P), pres- 
sure gradient, pore velocity (u), Reynolds number (Re), heat 
transfer coefficient (h), and surface area per unit volume (2). 


As previously stated in the model development, these system 


properties are functions of temperature, and hence, change 
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with time. In addition, since surface recession may occur, 
the spatial thickness (L) of the porous medium is provided 


as output data. 


D. EFFECT OF GEOMETRIC PARAMETERS 

In the initial effort by Vatikiotis [2], it was observed 
that, for a given set of boundary and initial conditions, the 
internal geometry of a porous medium significantly influenced 
combustion. Results were reported describing the dependence 
of combustion on particle size, and the thickness of the 
porous medium. As a refinement of that initial work, results 
are presented below which show the influence of the geometric 
parameters (i.e., permeability, porosity, porous medium thick- 
ness) on the minimum initial temperature needed to sustain 
combustion. "Minimum initial temperature" is defined here as 
the minimum constant-value temperature distribution (i.e., 
slope equal to zero) that will sustain combustion when used 
as an initial condition. A lower temperature will cause the 
porous medium to cool to ambient temperature. As will be 
Shown in Section V.E, both the shape and magnitude of the 
temperature initial conditions influences the combustion proc- 
ess. In order to form a basis for comparison, the analysis 
was performed with uniform distributions for the initial con- 
ditions (the carbon and air having the same initial tempera- 
ture). The initial condition for the oxygen concentration 
(also uniform) was taken as the concentration found in air at 
the initial condition temperature and ambient pressure. However, 


any Oxygen concentration distribution, as an initial condition, 


79 





would not have changed the results. The boundary conditions 
used in the analyses were the insulated boundaries, expressions 
TII.45 and III.46, for the carbon temperature, and the Danck- 
werts conditions, expressions III.47 through III.50, for the 
air temperature and oxygen concentration. The results of the 
analysis were obtained by varying the initial conditions until 
the "minimum initial temperature" for sustaining combustion 
was obtained. A geometric parameter was then changed to a new 
value while all other parameters remained the same. A "mini- 
mum initial temperature" for this problem was obtained. This 
procedure was repeated for several values of the same geometric 
parameter while keeping all other parameters fixed. The re- 
sults obtained showed the "minimum initial temperature" which 
sustained combustion as a function of the respective parameter. 
The results for each geometric parameter are discussed 
individually. 
ini tects Of Permeability 

Figure V.9 shows the dependence of the "minimum initial 
temperature" (defined previously in this section) on permea- 
bility. The description of the porous medium and the ambient 
conditions used in the following analysis is given in Table 
V.3. The range of permeability was determined by the magni- 
tude of the Reynolds number. The largest permeability 


(m = 4x 107° 


eee). at a pressure differential of 1.5 p.s.l1., 
resulted in an average Reynolds number through the porous 


medium of approximately 5. It was felt that extending the 
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TABLE V.3 --- Geometry of porous medium and ambient conditions 
for the permeability analysis in Section V.D. 


Particle Shape: spherical 
Particle Diameter (d): (various) 

Unit Cell Thickness (0): (various) 
Spatial Thickness of Porous Medium (LL): 1.0 in 

Porosity (Pp) » 476 
Permeability ‘M™) (various) 

Bulk Thermal Conductivity of Carbon (K,): 86.0 Btu/ft-hr-F 
Bulk Specific Heat of Carbon (GG): peor OG, bina e 
Bulk Density of Carbon ((,): Ula & Lbm/fto 
Thermal Emissivity of Particles (€): oy 

Ambient Temperature (T,): 80.0 deg-F 
Ambient Pressure (P): 14.7 psi 
Pressure Differential (AP): (various) 
Ambient Oxygen Concentration (Q,): 0172 lbm/st> 
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analysis to larger values of Reynolds number would exceed the 
limitations of the combustion model. Specifically, the 
assumption that Darcy's law governs the internal flow for 
Reynolds numbers greater than 5 would be questionable. Simi- 


9 ft‘), ata 


larly, the smallest permeability (m= Slo2 x on 
pressure differential of .1 p.s.i., gave an average Reynolds 
number of approximately .005. The lower limit Reynolds num- 
ber for which Darcy's law remains valid is not known. There 
is little information in the literature that addresses a lower 
limit Reynolds number. Although lower Reynolds numbers may 
be valid, the analysis was restricted to a minimum Reynolds 
number of .005. As can be seen in Figure V.9, the "minimum 
initial temperature" for sustaining combustion 1s a monotonically 
increasing function of decreasing slope of permeability. In 
other words, as permeability increases, the temperature needed 
for sustained combustion also increases. Important to the re- 
sults are that pore velocity also increases with increasing 
permeability. As stated previously, insulated boundaries 
were imposed on the porous solid. As will be shown in Section 
V.E, the results change when heat transfer occurs at the sur- 
faces of the porous solid. 

ao Strects of Porosity 

Figure V.10 shows the dependence of the "minimum 

initial temperature" for sustaining combustion on porosity. 
The description of the porous medium and the ambient conditions 


used in the following analysis is shown in Table V.4. The 
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TABLE V.4 ~~~ Geometry of porous medium and ambient conditions 
for the porosity analysis in Section V.D. 


Particle Shape: spherical 
Particle Diameter (d): (various) 

Unit Cell Thickness (0D): (various) 
Spatial Thickness of Porous Medium (L): Pe )e ote 

Porosity (Pp) (various) 
Permeability (mM) (various) 

Bulk Thermal Conductivity of Carbon (k,): 2G Os iehe yi get SI gaat as 
Bulk Specific Heat of Carbon (C,.): 225.  Seu7l om-F 
Bulk Density of Carbon (Pd 2 Weg s&s ibm/tt™ 
Thermal Emissivity of Particles (€): Sed 

Ambient Temperature (|): 80.0 deg-F 
Ambient Pressure (PP): 14.7 psi 
Pressure Differential (AP): =a oot 

Ambient Oxygen Concentration (Dos AOI gee lbm/ft? 
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results show that as porosity increases, the temperature 
needed for sustaining combustion decreases monotonically. 
This behavior is opposite, but less pronounced than that ob- 
served for the permeability. In addition, pore velocity de- 
creases as porosity increases. It seems reasonable that an 
increase in porosity provides more oxygen per volume of por- 
ous medium for combustion. The lowest value of porosity in 
the analysis (p = .48) was restricted by the idealized geometry 
shown in Figure III.1. For the idealized geometry, the low- 
est porosity occurs for the ratio of particle diameter to 
unit cell thickness, ad/D, equal to 1. The maximum value of 
porosity in the analysis was limited by the highest porosity 
reported for spherical particles in a stable configuration. 
As stated by Scheidegger [9], this value is .875. 
3. Effects of Porous Medium Thickness 

Figure V.1l1 shows the dependence of the "minimum 
initial temperature" for sustaining combustion on the porous 
medium thickness. The description of the porous medium and 
ambient conditions used in the analysis is shown in Table 
V.5. The range of thicknesses (.25"-4.0") investigated was 
limited (as in the case of permeability) by the Reynolds num- 
ber. As seen in Figure V.1l, the "minimum initial temperature" 
1s a monotonically decreasing function with decreasing nega- 
tive slope of porous medium thickness. In other words, as 
the porous medium thickness decreases, the temperature needed 


to sustain combustion increases. In addition, the pore velocity 
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TABLE V.5 --- Geometry of porous medium and ambient conditions 
for the thickness analysis in Section V.D. 


Particle Shape: spherical 
Particle Diameter (d): (various) 

Unit Cell Thickness (QD): (various) 
Spatial Thickness of Porous Medium (L): (various) 
Porosity (Pp) ~ 476 
Permeability (mM) (various) 

Bulk Thermal Conductivity of Carbon (k,): 86.0 Btu/ft-hr-F 
Bulk Specific Heat of Carbon (C,): foo) Seu? Lom - 
Bulk Density of Carbon (Dd: Toes lbm/ft> 
Thermal Emissivity of Particles (€): a] 

Ambient Temperature ( IT): 80.0 deg-F 
Ambient Pressure (P)): ae] Dis 1 
Pressure Differential (AP): —-.o psi 

Ambient Oxygen Concentration (D) = OL 72 lbm/#t> 
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decreases as the thickness increases. As in the case of the 
permeability, the results showed a wide range of temperatures. 
The "Minimum initial temperature" for sustaining combustion 
may vary as much as 700 degrees Fahrenheit over a range of 


thicknesses for a particular porous medium. 


E. EFFECTS OF EXTERNAL PARAMETERS 

A similar analysis, as discussed in the previous section, 
was performed for the pressure differential. The boundary 
conditions for the carbon particle temperature equation were 
then changed from the insulated boundaries, expressions III.45 
and III.46, to the radiation heat transfer boundary conditions, 
expressions III.5l1 and III.52. This permitted heat transfer 
from the surfaces of the porous solid. The emissivity, €, was 
assumed equal to .9. With the new boundary conditions, the 
effects of permeability and pressure differential on the "mini- 
mum initial temperature" for sustaining combustion were again 
analyzed. 

1. Effects of Pressure Differential 

Figure V.12 shows the effects of pressure differential 

on the "minimum initial temperature" for sustaining combustion. 
The description of the porous medium and the ambient condi- 
tions used in the following analysis is shown in Table V.6. 
The range of pressure differentials considered (.1-2.0 p.s.i.) 
was governed by the Reynolds number for which Darcy's law re- 
mains applicable. This point was discussed in Section V.D. 


In the initial work of Vatikiotis [2], the maximum pressure 
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TABLE V.& --- Geometry of porous medium and ambient conditions 
for the pressure differential analysis in 
Section V.E. 


Particle Shape: spherical 
Particle Diameter (d): (various) 

Unit Cell Thickness (0): (various) 
Spatial Thickness of Porous Medium (L): PIOeEH 

Porosity (P ) »~476 
Permeability (™) (various) 

Bulk Thermal Conductivity of Carbon (k,): 86.0 Btu/ft—-hr—-F 
Bulk Specific Heat of Carbon (C,): e201 Btu/1lbm—-F 
Bulk Density of Carbon (P.>: 19535 ieee 
Thermal Emissivity of Particles (€): 19 

Ambient Temperature (T)): 80.0 deg-F 
Ambient Pressure (P): 14.7 psi 
Fressure Differential (AP): (various) 
Ambient Oxygen Concentration (PD): nol lbm/ft> 
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differential considered was approximately .35 p.s.i. As seen 
in Figure V.12, the "minimum initial temperature" is a mono- 
tonically increasing function of pressure differential. The 
results show that, depending on pressure differential, the 
Minimum temperature for sustaining combustion could vary as 
much as 400 degrees Fahrenheit for a particular porous medium. 
In addition, pore velocity increases as pressure differential 
increases. 
2. Effects of Boundary Conditions 

In this section, the results of allowing heat transfer 
at the boundaries of the porous solid will be compared to the 
results obtained with insulated boundaries. Figures V.13 and 
V.14 show the results of the permeability and the pressure 
differential analyses, with and without heat transfer from 
the surface of the porous solid. The description of the por- 
ous medium and ambient conditions for the permeability and 
pressure differential analyses are shown in Tables V.3 and V.6, 
respectively. As stated previously, the radiation heat trans- 
fer boundary conditions were those shown by expressions III.5l 
and III.52. As seen in Figures V.13 and V.14, heat transfer 
from the boundaries significantly affects the "minimum initial 
temperature" at the lower values of permeability and pressure 
differential (or at lower pore velocities). The results show 
that when heat transfer occurs from the boundaries, the "mini- 
mum initial temperature" is no longer a monotonic function of 


permeability and pressure differential. This behavior was 
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suggested in the intial work of Varikiotis [2], where in addi- 
tion to radiation heat transfer, convection heat transfer 
was considered at the boundaries. The additional convection 
heat transfer causes the effects produced by permitting radia- 
tion heat transfer at the boundaries to become more pronounced. 
Moreover, these results are consistent with those reported 
Seevulis [13], Fontenot [51], and Spalding [52]. The reason 
for this behavior was summarized in Section V.C and will be 
explained in Section V.G. 
See ELnects of Initial Conditions 

The previous analyses of determining the effects of 
System parameters and boundary conditions were performed using 
constant-value initial conditions (1i.e., slope of temperature 
distribution was equal to zero). The constant-value tempera- 
ture defined the "minimum initial temperature" and was used 
as a basis for comparing the results. To also show the effects 
of the initial conditions on combustion, a single porous medium 
was subjected to numerous initial conditions. The analysis 
was performed as follows. The initial conditions of the car- 
bon and air temperatures were given a constant slope as shown 
in Figure V.15. With the slope fixed, the temperature distri- 
bution was shifted (1.e., moved up or down) until the thres- 
hold for sustained combustion was determined. That is, shift- 
ing the initial temperature distribution slightly lower would 
result in the porous medium cooling to ambient temperature. 
This procedure was repeated for several initial conditions of 


temperature having different constant slopes (both negative 
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FIGURE V.15 --- Constant-slope initial condition. 
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and positive). The boundary conditions used in the analysis 
were the Danckwerts conditions, expressions III.47 through 
TiIzr.50, for the air temperature and oxygen concentration, and 
insulated boundaries, expressions III.45 and III.46, for the 
carbon temperature. The results of the analysis are shown in 
Figure V.16. In Figure V.16, the average temperature is plotted 
as a function of AT defined by T(L) - T(0) for the particular 
temperature distribution. The average temperature, Tavg’ 1s 
defined as (T(L) + T(0))/2. As seen in Figure V.16, the curve 
is not symmetrical about AT = 0. This result is atrributed to 
the direction of the pore velocity. The largest difference 

in temperature between the porous solid and the air (and great- 
est heat transfer) is at the air inlet region (1i.e., x = 0). 

At the alr exit region (1.e., x = L), the temperature differ- 
ence between the porous solid and the air is small. The com- 
bined effect of high heat transfer at the air inlet region 

and low heat transfer at the air exit region produces the 
following result. There is a greater likelihood that sustained 
combustion will occur for a positive slope temperature dis- 
tribution than for a negative slope temperature distribution 
having the same average temperature and AT as defined above. 

In other words, for a positive slope initial temperature dis- 


tribution, sustained combustion occurs at a higher temperature. 


oer r ECTS OF PORE VELOCITY 
The way the parameters (i.e., permeability, porosity, 


pressure differential, and porous medium thickness) affected 
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combustion is in proportion to their effect on pore velocity. 
In other words, the results of the previous sections showed 
that pore velocity was affected by the parameters in the same 
Manner as the “minimum initial temperature" for sustained com- 
bustion (1.e., minimum temperature increased as pore velocity 
increased, and decreased as pore velocity decreased). An 
exception to this occurred when heat transfer was permitted 
from the boundaries. This suggests that pore velocity is the 
dominant variable affecting the combustion process. Since 
pore velocity is a function of the parameters analyzed (1l.e., 
permeability, porosity, pressure differential, porous medium 
thickness), it seems reasonable that the "minimum initial 
temperature" for sustained combustion may be posed as a func- 
tion of pore velocity. With this approach, the "minimum initial 
temperature” for sustained combustion of a porous medium with 
unknown properties could be determined by specifying the 
magnitude of the pore velocity. Moreover, pore velocity can 
be represented by specific volumetric flowrate of filter velocity 
(i.e., V/A) which is easily measured. Figure V.17 shows the 
"Minimum initial temperature", defined in Section V.D, as a 
function of the filter velocity at ambient temperature. The 
graph was made from the data obtained in the analyses of Sec- 
tions V.D and V.E (for insulated boundaries on the porous 
solid). The data in Figure V.18 has a spread of approximately 
+100 degrees Fahrenheit. The equation for the curve shown 


Per igure V.17 is, 
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where T. represents the "minimum initial temperature" dis- 


t 
tribution in degrees Fahrenheit, and Q is the filter velocity 
in ft/hr at an ambient temperature of 80 degrees Fahrenheit. 
Figure V.17 shows promise as a simple method for determining 
the "minimum initial temperature" for which sustained combus- 
tion can be expected in a porous medium. Perhaps by refining 


the formulation for calculating the pore velocity, there would 


be a smaller spread in the data. 


G. ANALYSIS OF RESULTS 

The Semenov model, discussed in Section III.C and illus- 
trated in Figure III.2, may be used to explain the results 
Mmemercea in Sections V.8B through V.F. Before continuing, cer- 
tain aspects of the Semenov model must be emphasized. First, 
the Semenov model does not govern combustion. It is an analy- 
tical tool which is used to describe the behavior of combustion. 
Secondly, the Semenov model as represented in Figure III.2 
is for a point in the porous medium. A "Semenov surface" 
would better explain combustion in the porous medium. This 
will be illustrated below. It follows from the Semenov theory, 
as discussed by Vulis [13], that "ignition" conditions (li.e., 
Sustained combustion) are determined by all the conditions of 
the combustion problem in a system. This is to say, that any 
change in system parameters (e.g., pressure differential, 


initial or boundary conditions) for a given porous medium 
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would result ina different Semenov surface. Lastly, the 
Semenov model is based on quasi-steady results. For example, 
the convection heat transfer, Tor and the heat generation, 
ea will be equal at their intersection (see Figure III.2). 
Conversely for the combustion model, at a particular point 
and time, Gd, may be smaller or larger than Bos depending on 
whether the temperatures are increaSing or decreaSing. How- 
ever, if the rates of change in the temperatures are small, 
the Semenov model can be used to analyze the results of the 
combustion model. 

Figure V.18 illustrates a "Semenov surface" generated 
from the results of the last example in Section V.C (Figures 
fey and V.8). In Figure V.18, Bs 1s a heat generation sur- 
face, and represents the particle temperature-heat generation 
relationship during the problem. The qo surface represents 
the convection heat transfer as a function of particle and 
air temperature (expression III.14). The particular heat 
transfer surface in Figure V.18 represents the instant the 
point at x/L = 0 transitioned from the kinetic regime to 
the diffusion regime of combustion. This can be seen by the 
tangent point, I, between the heat generation surface and the 
heat transfer surface. A difference between "classical" 
Semenov S-curves (Figure III.2) and those shown in the heat 
generation surface of Figure V.18 is the sudden drop in heat 
generation once a certain temperature has been reached. This 


1s only observed in the region away from the air inlet (x/L = 0) 
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FIGURE V.18 --- Semenov surface. 
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and is attributed to the total depletion of oxygen in that 
region of the porous medium. One might say this is a weak- 
ness in the Semenov model when applied to porous media. 
Conversely, it might be argued that when the oxygen is totally 
depleted the combustion problem becomes a heat transfer prob- 
lem (in a pure sense) of which the Semenov model has no 
application. Also in the region where the oxygen has been 
depleted, it 1s not necessary for dG, to intersect Ro: 
In order to visualize the results of the combustion model 
more clearly, the "Semenov surface" is abandoned for Semenov 
curves (illustrated in Figure III.2). Figures V.19, V.20, 
and V.21 show that Semenov model representations for locations 
x/L = 0, x/L = .5, and x/L = 1, respectively, of example C in 
Section V.C (transient results shown in Figures V.7 and V.8). 
The heat generation curves are shown along with heat transfer 
lines representing the film cooling at two different points 
mietime. During this discussion, the term “film cooling” will 
be used to mean the heat transfer from the particle to the 
air. "Convection heat transfer" will refer to the transport 
of energy by the internal flow. In addition to film cooling, 
heat transfer by conduction and radiation occurs within the 
porous medium. This has an effect of slightly curving the 
heat transfer line upward (1i.e., heat transfer becomes a non- 
linear vice a linear function of temperature). However in 
this problem, calculations showed that the ratio of film 


cooling to conduction and radiation was approximately 5, and 
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therefore, conduction and radiation will not be considered 
in order to simplify the graphical analysis. 

Comparing Figures V.19-V.21 shows that the maximum heat 
generation rates are different depending on the location in 
the porous medium. It was mentioned in the general discussion 
of the Semenov model (Section III.C) that implicit in the 
shapes of the S-curves are the effects of oxygen concentra- 
tion. In other words, decreasing the oxygen supply will 
cause the S-curves to peak at lower temperatures. Since it is 
more difficult for the oxygen to penetrate as the depth of 
the porous medium increase (1.e., oxygen is consumed as the 
flow passes through the porous medium), it is a reasonable 
outcome that the maximum heat generation rates decrease as 
the depth increases. The difficulty for the oxygen to pene- 
trate as the depth increases produces another effect as seen 
in Figure V.21. Unlike those in Figures V.19 and V.20, the 
S-curve and heat transfer curve does not have a tangent point 
which defines the “critical ignition condition". This be- 
havior is called noncritical combustion, and as discussed by 
Vulis [13], is characteristic of oxygen poor combustion. 
Moreover, these observations are supported by the experimental 
results of Thomas, Stevenson, and Evans [15]. In their 
experiments, it was shown by decreasing the ambient partial 
pressure of oxygen, that a critical combustion process (1.e., 
one in which a "temperature jump" occurs) transforms to a 


noncritical combustion process. 
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As seen in Figure V.7, the temperature first rises ina 
local region of the porous medium (in this case near x/L = 1) 
before moving aS a wave toward the air inlet surface. The 
direction of movement was a characteristic of all the prob- 
lems in which sustained combustion was observed. This behavior 
is explained as follows. Temperature rises in a region where 
the heat generation 1S greater than the heat transfer between 
the particles and the air. The region of higher temperature 
widens as heat is transferred by conduction and radiation. 

In addition, the increasing temperature simultaneously in- 
creases the reaction rate. At some moment, the oxygen enter- 
ing through the upstream side of the high temperature region 
is totally consumed within that region. When this occurs, 
the simultaneous increase in temperature and reaction rate 
wiil occur only on the upstream side of the region of high 
temperature. As the increasing temperature moves towards the 
air inlet surface, the front of the region of depleted oxygen 
also moves nearer the air inlet surface. This results in the 
temperature response appearing to move in the form of a wave 
toward the air inlet surface. It is interesting to note, that 
the magnitude of heat transfer by conduction and radiation in 
the increasing temperature region waS approximately twice 
that of the convection heat transfer at the outset of the 
example (Figures V.7 and V.8). It seems reasonable that by 
increasing the pore velocity to a sufficient magnitude, the 


region of increasing temperature would not propagate but blow 
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emt (L.€., toward the x/L = 1 surface). In fact, this is 
observed for the problems in which the porous medium cools 
to ambient temperature. 

Up to this point, the thermal response of the porous medium 
has been discussed in terms of "Sustained combustion" or 
“cooling to ambient temperature". The Semenov curves gener- 
ated from the results of the combustion model can be used to 
determine transition from kinetic to diffusion combustion. 

The kinetic and diffusion regimes of combustion are discussed 
in Section III.C. In a practical sense, the exothermic proc- 
ess occurring in the kinetic regime is synonymous with the 
term "oxidation". The process occurring in the diffusion 
regime is simply called "combustion". From everyday experi- 
ence, one knows that oxidation can be a slow process, occurring 
at low reaction rates. Conversely, combustion can produce 
considerable heat, reflected by high temperatures and reaction 
rates. As will be shown, transition from oxidation to combus- 
tion occurs rapidly. In addition, the porous medium may be 
undergoing oxidation in one part, and combustion in the other. 
It 1S important to understand the behavior of the thermal 
process at this level. Especially, if there 1s a need to 
closely control the process. 

Each point within the porous medium will transition from 
Oxidation to combustion at different times. This can be 
seen in Figures V.19-V.21 which were obtained from the results 


of example C. The locations x/L = l, x/L = .5, and x/L = 0 
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transitioned at 12.5 seconds, 22 seconds, and 41 seconds, 
respectively. Therefore at the start of the problem, the 
porous medium as a whole, was undergoing oxidation. More- 
over, transition from oxidation to combustion started at 

x/L = 1 and moved to x/L = 0 with time. This is consistent 
with the behavior observed for the temperature response. As 
discussed by Vulis [13], the theoretical value for "critical 
ignition" temperature of the carbon can be calculated using 


the "N. N. Semenov equation" given by, 


= aR tt : (1 = 4RT/E) 1/7) (v2) 
where E is the activation energy of the Arrhenius expression, 
Se is the universal gas constant, and 7 is the absolute air 
temperature. This equation is derived by equating the 
Arrhenius expression III.16 to the convection heat transfer 
expression III.14, and requiring the slopes be equal. Note 
that the theoretical "critical ignition" temperature is a 
function of air temperature and activation energy only. Since 
the transition from oxidation to combustion at x/L = 1 shows 

a noncritical behavior, equation V.1l1 does not apply. The 
theoretical "critical ignition" temperatures for locations 

x/L = 0 and x/L = .5 are 1831 degrees Fahrenheit and 1492 
degrees Fahrenheit, respectively. The values obtained by the 
combustion model (Figures V.19 and V.20) are approximately 
1825 degrees Fahrenheit for x/L = 0, and 1440 degrees Fahren- 


heit for x/L = .5. Once the “critical ignition" temperature 
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was reached, the increase to the "combustion" temperature 
occurred in less than 1 second for both cases. Examples A 

and B in Section V.C (Figures V.2-V.6) were the same problem as 
example C, discussed here, except combustion was initiated 
differently. Example C was initiated with the porous medium 

at a constant temperature known to sustain combustion. [In 
examples A and B, the porous medium was subjected to a heat 
flux at the air inlet surface (x/L = 0) with the porous medium 
at ambient temperature. The amount of time for which the heat 
flux was applied determined if sustained combustion occurred. 
Sustained combustion was observed for example A in Section 

V.C. The amount of time the heat flux was applied in example 

B was not sufficient for sustaining combustion and the porous 
medium cooled to ambient temperature. It is interesting to 
note that the temperature at x/L = 0 in example A reached the 
"Critical ignition" temperature of 1825 degrees Fahrenheit. 

The maximum temperature obtained at x/L = 0 in example B was 
approximately 1650 degrees Fahrenheit. There is good agreement 
between the theoretical values of "critical ignition" tempera- 
ture and those obtained by the combustion model. 

The pore velocity significantly affects the combustion 
process. This was discussed in Section V.F. Increasing pore 
velocity has two effects. First, the internal heat transfer 
coefficient, h, is increased, thereby, increasing the amount 
Of film cooling. Secondly, the supply of oxygen is made greater 
(i.e., convection mass transfer of oxygen molecules increases). 


The combined effect is that by increasing pore velocity, higher 
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local temperatures are needed for transition of the porous 
solid from oxidation to combustion. In other words, assume 
that transition occurs at a particular air temperature fora 
given pore velocity. Increasing the pore velocity without an 
increase in air temperature may not be sufficient for transi- 
mem tO occur. This is illustrated by Figure V.22. Increas- 
ing the pore velocity will increase the slope of the heat 
transfer line, and also change the shape of the heat genera- 
tion curve as shown. It iS apparent that an increase in air 
temperature is needed to make the steeper-sloped heat transfer 
line tangent to the heat generation curve. In addition, higher 
pore velocities produce higher combustion temperatures once 
transition occurs. In the general discussion of the Semenov 
model, it was pointed out that the temperature dominates the 
combustion process in the kinetic regime. Increasing the 
oxygen supply will have the greatest effect in the diffusion 
regime (1.e., increasing the combustion temperature). This 
behavior described by the Semenov model, combined with the 
effects produced by convection heat transfer (energy trans- 
port by internal flow) and to a lesser extent, heat transfer 
by conduction and radiation, will determine if the porous 
medium will undergo sustained combustion for an increase in 
pore velocity. 

Figure V.23 shows the effects of decreasing the pore 
velocity. The slope of the heat transfer line becomes smaller, 


and the heat generation curve changes shape as shown. [In 


iS 
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PIGURE V.22 --- Effects of increasing pore velocity. 
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FIGURE V.23 --- Effects of decreasing pore velocity. 
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particular, lower combustion temperatures will occur in the 
diffusion regime. Unlike the increasing pore velocity case 
above, there is less certainty in predicting the behavior of 
combustion when pore velocity is decreased. It appears that 
conduction and radiation heat transfer become more of a factor 
at the lower pore velocities. This seems reasonable since the 
ratio of conduction and radiation heat transfer to convection 
heat transfer increases with decreasing pore velocity. In 
addition, the effects of heat transfer from the boundaries at 
the lower pore velocities must also be considered (shown in 
Figure V.13 and V.14). The combined effects of all the heat 
transfer mechanisms, including the boundary conditions, will 
determine if the porous medium will undergo sustained combus- 
tion as the pore velocity is decreased. 

Lastly, the effects of the boundary conditions must be 
considered separately. The results of changing the boundary 
conditions were presented in Section V.E. It was observed 
that for heat transfer occurring at the surfaces of the porous 
solid (e.g., radiation boundary conditions), higher initial 
temperatures were needed to sustain combustion. Allowing heat 
transfer from the surfaces of the porous solid (either by 
convection or by radiation or by both), affects the combustion 
process in the boundary condition regions as follows. Return- 
ing to the Semenov model, there is nothing explicitly asso- 
Clated with the boundary conditions that would change the shape 


of the heat generation curve. In addition, since pore velocity 


L6G 





determines the magnitude of the internal heat transfer coeffi- 
cient, it appears that the slope of the heat transfer line 

will also be unaffected. However up to this point, only film 
cooling was considered in the heat transfer line of the Semenov 
model. Noting that heat transfer from the porous solid boun- 
daries results in steeper temperature gradients in regions 

near the boundaries, heat transfer by conduction and radiation 
must also be considered. With the additional components of 
heat transfer, the heat transfer line will curve upward. MThis 
means that the heat transfer changes from a linear function 

to a nonlinear function of temperature. The change of the 
shape is illustrated in Figure V.24. As can be seen, the 
increased heat transfer by conduction and radiation may pre- 
vent transition from oxidation to combustion. This, in turn, 
will restrict the heat generation to the kinetic regime. It 

is reasonable to expect that as the reaction moves towards the 
air inlet boundary, the combined effects of greater heat trans- 


fer and lower heat generation will oppose sustained combustion. 
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FIGURE V.24 --- Effects of boundary conditions on 


heat transfer. 
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Vi CONCLUSIONS 


The main objective of developing a model which predicts 
the combustion behavior of a porous medium has been achieved. 
The results of the mathematical model are in good agreement 
with those obtained by experimental methods. Moreover, the 
analyses presented have shown that combustion and heat trans- 
fer in porous media is a complex process involving the inter- 
action of heat transfer and heat generation. This interaction 
depends upon the geometric parameters of a porous medium, 
boundary conditions, environmental and initial conditions. 
The behavior of a system can change radically by altering any 
number of parameters. This point was demonstrated by examples 
A and B in Section V.C where a difference of one second in 
applying a surface heat flux meant combustion or extinguish- 
ment. The results have also shown that conduction and radia- 
tion between particles may play a significant role in deter- 
mining combustion behavior and should be accounted for. The 
computer program makes it possible to look at a large number 
of cases. Similar analyses by experimental methods would be 
economically impractical. 

Based on literature surveys performed during this inves- 
tigation, it appears that much of the engineering develop- 


ment associated with porous media involves a trial and error 
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process. It seems reasonable that an analytical model would 
be an essential tool for the engineer to employ in the de- 
Sign process. In this investigation, the model has been 
assessed for its applicability, and to some extent, its 
accuracy. It is hoped that in the future the model will be 
used to determine the effects of the design variables (l.e., 
the geometric parameters, pressure differential, etc.) on 
performance. Specifically, the model shows promise for 
evaluating the combustion efficiency and stability of a sys- 
tem. It is in this capacity that the combustion and heat 


transfer model will be of greatest value. 
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APPENDIX A 


FORMULATION OF FIELD EQUATIONS 


1. PRESSURE DISTRIBUTION EQUATION 


Darcy's law for one dimensional flow, neglecting body 


forces, is, 


_ mn aP 


vy Cae 


CAs Bb) 


Substituting in the Dupuit-Forcheimer relation, and solving 


for u, equation A.1l becomes, 


vy = - BGP) 
pu dx 


The continuity equation (derived in Appendix B) is, 


d(pp_) 
a 9 = 
me ox eae = 


Substituting equation A.2 into equation A.3 yields, 


(pp .) re) ae dap 
at ox Ur ax 


Z dO 
eS | oe eg 
dx Pa 


IPAs 


(A.2) 


(A. 3) 


(As) 


(A.5) 





Equation A.5 with associated boundary conditions (presented in 
Section II.B) is solved for the pressure and pressure gradient 
distribution. The pressure gradient is then substituted into 


equation A.2 to obtain the pore velocity. 


2. POROUS SOLID HEAT TRANSFER EQUATION 

To perform energy balances on both the porous solid and 
on the air, a differential volume of porous medium was segre- 
gated into respective volumes of constituents, that is, 
qv. = (l-p)dV for the solid, and dv. = pdV for the air (shown 
am Figure A.1). The convention used for the energy balance 


of an arbitrary differential volume, dV, is, 


Heat into Heat S Heart ouc it Increase in 
dv Generation of dav Internal Energy 


The heat transfer mechanisms considered for the carbon 
particles are conduction, radiation heat transfer between 
particles, convection heat transfer from the particles to 
the air, and heat generation. Applying the above convention, 
the energy balance on a differential volume of porous solid 
ZS, 


(1°P) deonaSA| + (l-p)q_44dA . ~ Ggen tA’ ENB ey) 


= (l-p)q dA | + (1-p)q_ da 


cond 
= X+dax X+dx 


+ | 
Slooeae = 


+ (l-p)q,_,dv 
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(t-9)-dV 
dv 2 


FIGURE A.1 --- Separating a differential volume of 
porous medium into respective volumes 
Oe SO lide andeead i, 
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Representing terms on the right side of expression A.6 by 
Taylor series expansions (neglecting higher order terms), 


the energy balance becomes, 


(1-P) deongtA| | + (1-p)q,,qd4 b a flog (A.7) 


dA 


= _ ore _ 
TS P) dona | a 7x b C1 P) doong! AxdA + (1 P) dag 


re) ; ; 
+ sy (l-p)a gl dxdA + qo, ydh' + (l-p) ai, av 


Subtracting terms, and rearranging, expression A.7 becomes, 


u ; 
gldV - S-[(1-p)q_,g]dv - qa 


Pe) 
mega! (-P)q conv 


con 


+ q naa ' = (1-p)q, ,dv (A.8) 


ge 


Substituting the following expressions into equation A.8, 


oT 
na #2 Ke => Fourier's law (A.9) 
oT A 
q = = 8 See Radiation analogy to (A.10) 
Back BES Fourier's law 
ony A(T. ~ T,) Newton's law (AI) 
Ioen = Ne Heat generation CATE eZ) 
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c— Internal energy AGS) 


yields, 
5 ae 
a —- eee ss a aA | t 
axel fl Pp) (k +k, ) ner J dv h(t. AN) GIR) RGA (A.14) 
aT. 
2 SEND Sa aay CU 


Dividing through by dV, and defining dA'/dV as z, the specific 
internal area (1.e., surface area per unit volume), equation 


A.14 becomes, 


oT 
at (1- emer 2 = 
mee. P) (k,+k_) aa hz(T. T,) + Be (A.15) 
: aT. 
= NED) Og Sey aes 


The expressions used to obtain the values of the properties 
and parameters in equation A.15 are presented in Section 


ia. E£. 


3. AIR HEAT TRANSFER EQUATION 

The formulation of the air heat transfer equation will 
begin with the general one dimensional energy equation, 

D 2 3 oe 


— a = 
Po, pele tau) = aylpki, a) + he(T -T,) + uF (A.16) 





As in the porous solid heat transfer equation, the time and 
position dependent porosity appears inside the differential. 
Expanding terms and neglecting body forces, the air energy 


equation becomes, 


ee 


D(>=u ) - oT. 
Meeipe 8 Pa peCO:C«COe Pe Kan ye? F BAIT -T,) et) 
= ax (PUP) = 7, (Put) 


Consider the momentum equation for the x-direction (neglecting 


body forces), 


3 = 2 o ea oe 
xp (PP,u) = 7, (PP, U) Tra ares, 7x (PP) (A.18) 
and the continuity equation, 

3 ee 5 One: = gu 

nes Onl & Ox PP! PP. ox (A.19) 


Multiplying the continuity equation through by u, and substi- 
tuting this into equation A.18, the momentum equation becomes, 
du du 3 


= 2 ony 2 So 
2 Ain ROE ae PE) — ae tPt.) (A.20) 


Q? 


Multiplying equation A.20 by u and noting that, 


du du 
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equation A.20 becomes, 


Du e) 0 
up P. ne = = 5x (PP) =U ee pines) Chee 2) 
The energy equation A.17, after substituting 
D(su*) /Dt = uDu/Dt and expanding terms, is, 


aT. 
aaa) + hz(T_-T) CN 2 3) 


Substituting equation A.22 into the above energy equation and 


cancelling terms, equation A.23 becomes, 


Es oY 0 a 5 OE or 
PPabt 9x ‘P Xa DX cca PE ox) Pf xx ox : 


The viscous dissipation term, PT, dU/ 9X, is neglected because 
the fluid is a gas flowing at a low velocity. Therefore, the 


energy equation for the air in the porous medium is, 


aa 


Pe... a _ du 
mee 9) 8x “a Gx? * DAITQT,) ~ PP aE eee 


With specific enthalpy for a gas defined by, 
h =e+ P/p (A.26) 
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pe/pt can be expressed as, 


pe _ Dh_ 1 pp, P a (A. 27) 
Dt Dt oy Dhe 2 te 
Ps 


a 


Multiplying the continuity equation through by P/ (po *) gives, 


Dp 
E a Pi De jg) el El 
0° pt po, Dt 0, 9x 


Substitution of the above expression into expression A.27, 
the substantial derivative of internal energy can be expressed 


as, 
ee ee ee (A. 29) 


Substituting expression A.29 into equation A.25, the energy 


equation reduces to, 


aT 
cigs , DP _ ee _ 2 a 
Meee - PO Pe = SQ (ke) + h2(t,-T,) (A. 30) 


The specific enthalpy can be represented in terms of tempera- 


ture, as, 
dh = Tds + - ap (Aes 1) 


where s is specific entropy. For a perfect gas, ds may be 


expressed by, 
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a 


ds = Sp ci oT 


dP (A 32) 


Substituting expresSion A.32 into expression A.3l1, and can- 


celling terms, enthalpy can now be written as, 


dq = cy aT (Re 33) 
Or, 
DT 
Dh _ a 
EW ca: SDE (A. 34) 


Substituting expression A.34 into equation A.30, the air 


energy equation or heat transfer equation becomes, 


5 dT. dT. 
ax (P Kk. “a = A See + HZ ae) (A. 35) 
dT 
D ve a 
eee 2 aca Ot 


Initial results showed that pP changed slowly with time. 


Therefore, the substantial derivative of pP, as shown here, 


D _ 2a 3 
De (PP) = aE (PP) a U 7x (PP) (A. 36) 


is reduced to ud(pP)Ax. The expressions used to obtain the 


properties and parameters in the coefficients of equation 


A.35 are presented in Section III.E. 
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4. OXYGEN MOLECULE DIFFUSION EQUATION 

The final consideration in formulating the field equa- 
tions for the model is the transport of oxygen molecules. 
The oxygen molecule transport equation is obtained by a con- 
servation of species belance on the differential volume of 
air, dv. = pdV. The convention used for the species balance 


is given by, 


0. mG _ Of out O05 . n O04 | 
dv of dv Consumption Accumulation 

The transport mechanisms considered were diffusion resulting 

from concentration gradients, convection, and consumption of 

oxygen by combustion. Applying the above convention, the 

species balance on the oxygen becomes, 


dAl +opm ga | = dA (A. 37) 


P™Maift 
da | +m dA' + pm, av 


cons CC 
xXx+dx 


Representing the terms on the right side by Taylor series 


expansions (neglecting higher order terms), the species balance 


becomes, 
PMaig¢¢ GA; + P Moony FA| PMs; ¢,¢ dA (A.38) 
xX x x 
- ° ° oo ; . ° 
: ox (PM ai ge) AXGA i PMeonyA| + 9% (PMcony) dda a Moons“ pe accay 
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Cancelling terms and rearranging, equation A.38 becomes, 


_ = fhe Bane 
eee aqice oY yx Pm 


Conv 


= ' = 
} av Monae PMocec 


(A.39) 


Substituting the following expressions into equation A.39, 


Mace = ~% = 
Nconv = U? 
“cons — EG 
2 
a ~ +e 
yields, 


2-(p 0, 22) av - Sup o)av - 


e gx ox 


Fick's law 
Convection transport 


Consumption by 
combustion 


Accumulation 


RecA 


Z ao 
oe), 1S) es oo) 


(A.40) 


(A.41) 


(A.42) 


(A. 43) 


(A.44) 


Dividing both sides by dv, and letting dA'/dV equal the 


Specific internal area, z, the oxygen molecule diffusion 


equation becomes, 


o 
ax ‘P Ue ax 


The methods and expressions for obtaining the properties 


De, - 
mea) 855? 


ll 
O 
ate 


(AAS) 


and 


parameters in the coefficients of equation A.45 are presented 


fm Section II.F. 
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5. TRANSFORMATION OF FIELD EQUATIONS FROM A FIXED COORDINATE 
TO A MOVING COORDINATE SYSTEM 


As discussed previously, the reaction rate at the air in- 
let surface of the porous medium results in surface recession. 
To account for this surface recession, the field equations 
must be transformed from a fixed coordinate to a moving coor- 
dinate system. The method of transforming the field equations 
will be shown for only the porous solid heat transfer equa- 
tion since the approach for the other field equations (i.e., 
air heat transfer, and combined Darcy's law and continuity 
equation) is identical. 

The porous solid heat transfer equation with the modified 
heat generation term is, 


oT 
Pe) Cc rod = 


oT. 
= (CSO eae 


Since the x coordinate is a function of time during surface 
recession (e.g., T (x(t) ,t)), the time derivative term in 


equation A.46 must be expanded using the chain rule, 


oT oT 
e = c dx o 
1 Es wh x 


(A.47) 


The x coordinate in the field equations is nondimensionalized 


by n= x/L. Substituting n for x, expression A.47 becomes, 


iL 





(<8) (S%) + 8) = —S) (32) (S%) + (8) 
t i x i t 3 
Noting that, 
OM. a oe 
gx L 
and, 
ad a dn Cir dL 
ae ~ “ae * ae ~ "Gt 
expression A.48 becomes, 
Pc, ax, ve _ n (aby 27 oe 
ox, dt’. gt y Lidt’ on dt 


(A.48) 


(A.49) 


(P50) 


(Aso) 


Upon substituting n = x/L into the remainder of equation A.46, 


and substituting expression A.5l for the time derivative term, 


the porous solid heat transfer equation for the surface 


recession problem is, 


me et (1- ) (k_+k oS) het = ) RUB = Q) 
an P ee on eo Bel g ! 


oT. ai? 


= _ i SC, Sl 
= (l-ple,c lela + ze 


As stated in Section III.G, the thickness, L, 


of time, and L can be obtained from expression III.44. 


CRB) 


as a function 


The 


air temperature and combined Darcy's law and continuity 


i3e 





equation appear similar to equation A.52 for the surface 


recession problem. 
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APPENDIX B 


FORMULATION OF AUXILLIARY EQUATIONS 


lL. CONTINUITY EQUATION 

The law of conservation of mass requires that the substan- 
tial derivative of the fluid mass in a differential volume be 
zero. Therefore, the continuity equation for a fluid in a 


porous medium is expressed by, 
D(pp_av) = 0 (B.1) 
Dt a ? 


or in an equivalent form, 


3 e) du 
7" (PP 4) + u aa (PO) + Pile aa = 0 (B.2) 


2. RADIATION HEAT TRANSFER ANALOGY TO FOURIER'S LAW 
Radiation heat transfer in the model was represented by 


an analog to Fourier's law of conduction heat transfer shown 


Day ; 


SC 
Grad ~ “Ke Oe 2 ac)) 
where Kk. 1S an equivalent conductivity of the particles due 
to radiation. The development is as follows. Assuming air 
to be transparent to radiation, and treating the idealized 


geometry of the porous medium (Figure III.1) as a series of 


iD 





closely spaced walls, the net radiation heat flux between 


adjacent walls is, 


= eC sc ae 64 
Trad ~ rae we ae (B.4) 
where T_ and T are the respective absolute wall tempera- 
x x+dx 


tures, € is the emissivity of the carbon, and o is the Stefan- 
Boltzmann constant. Expanding cine in a Taylor series about 
Beegx' and neglecting higher order terms, the series expan- 


Sion may be written as, 


_ G E(t ‘ Tt) _ 4ce 7 dT 


drad 2——°°X x 2-— “xt+dx Ox dx (B.5) 


Simplifying, the above expression becomes, 
dx (B. 6) 


Equating expressions B.3 and B.6, and noting that 93T/dx = dT/dx, 


k_ becomes, 

ie 

(Bon) 
where dx is now equal to 6, the pore diameter. From the close 
Spacing of the particles, the temperature difference will be 


small as compared to the magnitude of the temperature. Noting 


this, the average absolute temperature of the carbon particle 
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may be substituted for T The equivalent radiation con- 


x+ax~” 


ductivity expression becomes, 





ie = 4oe6 m3 (B.8) 


and the radiation heat transfer from particle to particle 


may be represented by, 





Reo 
ee eee) 
q 7 Tae a aoe) 


For small pore diameters equation B.9 will be a good approxi- 


mation of the radiation heat transfer between particles. 


3. POLYNOMIAL APPROXIMATIONS OF THERMAL PROPERTIES 
Relations giving the dynamic viscosity, thermal conduc- 
tivity, and specific heat at constant pressure of air at 
different temperatures were required. A simple method to 
obtain values for these properties is to fit empirical data 
with 2nd order Lagrange polynomials. 
The general form of the 2nd order Lagrange plynomial is, 
(T,-T,) (T,-T,) (T,-T,) (T,-T,) 


Ki= Tpopypmepy *1. * tropa aas*~ CX (B.10) 
- (Ty To) (Ty T.) 1 (T. T,) (8, T,) 2 


(T,-T,) (T,-T,) 


ek 


where Ks 1s the property value at the ith temperature, T, 


ey 





Choosing three temperatures that are representative of 
those observed during the analysis, the corresponding values 


of the properties are, 


(subscript) Temp. u Ke Cc. 
deg-F lbm/ft-hr Btu/ft-hr-F Btu/lbm-F 

i So «(MN P873x10 - 1.5161x10 - 24020 

2 Revome i soo4stoe-  3-9013sx10 - 27268 

3 EL OMOIGEPG on = = ««7.1G¢47x110 mo 


Applying expression B.10 to each set of properties results 


in the following set of polynomials, 


u o= 2 208) % Loe 18 i 4.633x 10-7, + 4.427x10 ° iB) 

k, = BeGos x Lou.” i + 1.930 x 107°T_ MeesGa cro) — aiBeden 
: Por -5 

c, = -1.293x10°7 To + 2.758x10 "rT, + .238 (B.13) 


Each expression gives property values within two percent of 


the data for temperatures to 3000 degrees Fahrenheit. 


4, JUSTIFICATION OF THE DANCKWERT'S BOUNDARY CONDITIONS 
The Danckwerts' boundary conditions for the air heat 


transfer are, 


k —-<-— = oO uci (T, - T,) (B.14) 
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= se (3. 3) 


Bischoff [41] presents a discussion of these boundary condi- 
tions applied to mass diffusion equations in porous media. 
An analogous discussion is made here for the Danckwerts' 
boundary conditions applied to the air heat transfer equation. 
The porous medium with entrance and exit regions is shown 
in Figure B.l. The air temperature, To for each section of 
the region is distinguished by subscripts, as are the proper- 
ties. The air heat transfer equations for each of the regions 


are as follows, 








a°r aT 
ay ay 
Ke ¥ = P. BT aan 7s = 0 x = 0 (B.16) 
dx 
A aT. aT. 
ax (© *. sax) moc ae hz(T-T.) a0 0 <  X < -L (Bley) 
aT aT 
ae ay) 
kK. : ae pP, Us Cc. Sasa 0 x > L (Beis) 
Se 
with the boundary conditions, 
T. ee =) t. (B.19) 
1 
a) = T (0) (B.20) 


ales, 





Flow e> 










Entrance 


Region Porous Slab 


age l< 8 
ara 


FIGURE B.1 --- Justification of the Dankwerts 
boundary conditions. 
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d = eS 

Mees adx' a, °°) = po uc ,T. (9) Pkigx!T, 9) 1 
(3252)) 
et = T(t) (B.22) 

d = —_ 

Seee-a a,’ ~ “aax'Ta,‘b)! = 29 fe EY) PKaS—(T_ (L) | 
(B23) 
me {o) = finite (B.24) 


The above boundary conditions impose the restriction that 
there is no convection heat transfer from the porous solid 
to the air at the boundary surfaces. The difficulty of mixing 
Danckwerts’ and convection heat transfer boundary conditions 
is discussed in Section III.H. For convenience, the thermo- 
physical properties will be treated as constant in the entrance 
and exit regions. 

An analytical closed-form solution of this set of equations 
is unknown because of the nonlinearity of the temperature 
dependent properties in equation B.17. However, the solutions 


of equations B.16 and B.18 are, 


oe el ae: 


iS a a 
i, eee ha exp (—— x) <0 (Be25) 
a 
T Oe oO) el 
a, = K,+ kK, —— x) x > L (B.26) 
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Applying boundary conditions B.19 and B.24, the following 


results are obtained, 
K = T (B27) 
K = 0 (Baz 3)} 
The solutions, B.25 and B.26, become, 


oe al Ca 
ieee it, + Ko exp(—~— 
if a 


xX) (B.29) 


* 
[A 
oO 


rF 
il 
~ 
~ 
[Vv 
ES 


(B30) 


It would be necessary to have the solution for the nonlinear 
equation B.17 to solve for K. and K.. However, the constants 
need not be known to continue with the analysis. 


From equation B.29, 


ee. pecans K. (Bie sk} 
and, 
d _ Pavia 
= a = KS K, (B.32) 


Substituting these into equation B.21l gives, 
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_ da 
Me a a2 wParicat2 ~— PPatC,7, °°) Pk ax (Ta '0) | 


CB B89) 
Cancelling terms and noting u, = Up Caren, U,P, = U5P,) 
yields, 
d 
po, uci, = po, uc, T, (0) - Pisa agen. 60) ] (B.34) 
Rearranging, the Danckwerts' boundary condition at x = 0 is, 
k Aap Cn = p, uc (Tt, =) (B35) 


Noting that K, is a constant, the derivative of equation B.30 


s' 
Ss , 
d . 
Ag ts = 0 (B.36) 
2 
and, 
Cir (L)] = 0 (B.37) 
ax a. ‘ 


Sbustituting these expressions and equation B.22 into equa- 


prom 3.23, and noting u, = up 2e., U,P, = U5P,) gives, 


d 
po, uc, T(t) - Ke ay (T, (L) 1 = pp,uc, T. (D) (B o3.0)) 


Upon cancelling terms, the second Danckwerts' boundary 


condition becomes, 
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ct = 
ae) i 0 (2529) 


An important consideration for using the Danckwerts' 
boundary conditions is that it simplifies the analysis since 
equations Similar to B.17 may be solved independently without 
having to consider entrance and exit regions. 

5. PRESSURE DIFFERENTIAL CALCULATION RESULTING FROM EXTERNAL 

FLOW AT SURFACE X = L 

To obtain the pressure differential, AP, for simulating 
the conditions of Fontenot's experiments [1] (discussed in 
Section III.H), Bernoulli's equation was used. The following 
observations showed this to be a valid assumption. Schlichting 
mee) POints out that for the ratio of u/U_ in the range of 
.0001 to .01, the effects of "blowing" or "suction" on the 
potential flow over the external surface of the porous medium 
may be neglected. A typical value of u/U, for the model at 
which U, = 25 knots was .0028. For steady flow over a flat 
plate, the flow field outside the boundary layer may be 
described by Bernoulli's equation. This is a direct result 
of the Navier-Stokes equation. In addition, the pressure 
gradient across the boundary layer may be taken as zero. 
Therefore, Bernoulli's equation, 


y? 
+ = = constant (B.40) 


— 
0 | 
w |r 


may be used to obtain the pressure differential across the 


plate. The parameters in equation B.40 are defined as follows: 
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P, pressure; 0._, density of air; U, free stream velocity of 
the air over an external surface of the porous medium. For 
the model, the density of the air was approximated by the 


ideal gas law as, 


O = P/R. Sty (B.41) 
where ne is the gas constant of alr, ie is the absolute tem- 
perature of the air, and P is the pressure. Substituting 
equation B.4l into equation B.40 gives, 


Ra T3 i 
f dP + = = constant (B.42) 





Upon integrating, equation B.42 becomes, 


ie .. ae) ue = constant (B.43) 
or, 
; us ; us 
a Ta In(P,) fore = Re "a, in(P.) on (B.44) 
Letting, 
oe 2 0, fT =f =F, P, =P, U, = U, 


and substituting these into equation B.44, yields, 
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U, 
= 


R a IDS i Re an In(P,) oa (B.45) 
Rearranging, equation B.45 becomes, 
¥2 
Ro ue In(P./P,) = on (B.46) 


Taking the exponential of both sides, and solving for Pie 


results in, 


_ 12 é 
PL = ip eye | U,/2R. cme) (B.47) 


From the above expression, and noting that AP = ee ee AP 


may be expressed as, 
AP = Pl fexp(-U2/2R_T,) - 1] (B.48) 
Expression B.48 was used to approximate the pressure differ- 


ential across the porous medium when simulating the conditions 


of Fontenot's experiments [1]. 
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APPENDIX C 


NUMERICAL FORMULATION 


i. FINITE ELEMENT METHOD 

The solution of the system of nonlinear, coupled, partial 
differential equations given by equations III.32, III.38, and 
III.39, subject to boundary and initial conditions, was ob- 
tained by a Galerkin formulation of the finite element method. 
The solution of equation III.12 was obtained using a shooting 
method. 

a. Galerkin Formulation 

A Galerkin formulation of the Finite Element Method 

was used to obtain solutions of the porous solid and air 
energy equations, and the oxygen diffusion equation. A con- 
venient form of equations III.32, III.38, and III.39 was 


used in the formulation as shown by, 


~23 oT. o- 
iL DN i oN ee oe ale Se = (l-p) 0 cae 
(eae) 
aval oT 
a2 oO a =) a 7 
L an tP kK a PO. cul aap + h2(t. T.) (G22) 
aT 
=-ldQ _ a 
+ ub ame aos ea 
~2 9 9d) pnt a = = a9 

ie an Pe an L a7 (PU >) a Pe = P5E (ey) 
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where the spatial coordinate, x, was nondimensionalized by 
n = x/L. 

The closed domain (0,1) was partitioned into (n-1) 
contiguous elements of variable length Ley cael. sll 2 
This defines an n nodal point model. The three field varia- 


bles, Toe To @ were approximated by, 


tem,t) = Yiim,t) = ) Gin)s, (e) (E54) 

ft) = Yo(m,t) = ) Gin)s,(+) (G25) 

o(n,t) = Waln,t) = } Gin), (t) (E26) 
where Gs for i=l,...,n is a set of specified basis func- 
tions with local support, and the sets {9 785, 9 33 det. eee 
are the solution coefficients to be determined. The G were 
selected tosatisfy the condition Sep = 3 where the 


Kronecker delta, iy is defined by Si = 1 for i= j, and 


5% = 0 for iF j. Asa result, O40 O55 and 9 3 are the 


13 
values Var Vor 3 at the nodal points (i.e., a = Y,(n,,t)). 
Linear interpolation functions (shown in Figure C.1) 
were used as the basis functions. These are the lowest 
polynomial functions which provide the necessary function 
Somrinuity. 


As a measure of error, a residual function, ta, is 


defined for each field equation by, 
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1 2 5 4 n-2 n-l n 


FIGURE C.1 --- Linear shape functions used in the 
Galerkin formulation of the FEM, 
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r, (nt) = A, ty) Ts zl 


Il 


Uae Cer») 


where AG denotes the spatial operator of the ith equation. 


For convenience, the following convention for differentiation 


is adopted, 


ee eee (Ca:3)) 


2 (a a (C.9) 


“ae (C.10) 


Merertield equations C.1l, C.2, and C.3, the residuals are, 


gi 


Pe n 
L“((1-p) (ka+k,) J Gio, J'-hz J G, (8, -8,) (c.1L) 
2 in| i=l i “fi 


8 (Cei2) 


IES 





2 n = n 
= - t (p02 SLOqehye & (pu ) G,9, ) (Ger 3) 
1=1 zr 1=1 a 
taf _ 
- R. @Q Z G.9 = ic G.9 
Oo i ae as ees 


where the coefficients multiplying the response variables are 
themselves functions of the response variables, and thus, 

the equations are nonlinear. In accordance with the Galerkin 
method, the final system of ordinary differential equations 
was obtained by setting each residual, Las orthogonal to each 


basis function, G5, that is, 


i 1 tee orcas esp Tl 
f G.r,dn = 0 (Garay) 
0 J 3 = 1,2,3 


The 3n ordinary differential equations given by equations C.14 
retain the character of the original set of partial differ- 
ential equations. Thus, linear field operators transform to 
Matrix operators and nonlinear, coupled field operators become 
nonlinear, coupled algebraic operators. Incorporation of 

the boundary conditions resulted in 3n nonlinear coupled 


ordinary differential equations, 


A(t) ¥ (95,9578) te F(t) = B(t) wp (Ga 5) 


Subject to initial conditions, where season < sh) Mace lx, 


i2 


1k 





A is the operator associated with the field operator A, in 


~! 


expression C.7, and F is an excitation vector. 


Adopting the convention, 


nN 
<G,> {6} = ! Gace: (C.16) 


and applying the operation of expression C.14 with an integra- 


tion by parts on the second order derivatives gives, 


iL 
-2 ' iL 2 = | ] 
(G,}L “(1-p) (k,+k,)<G,> (8,31, -L “(1 P) (kath) f {G,} <G,>'dn 


1 1 
-hz ¥ (G;}<G,>dn{6,} + hz y {G;}<G,>dn{65] (eras?) 


i 


il iL : 
+ R263 % (Gj }<G,>dn{8,} = (1-p)o.¢, a (G;}<G,>dn{9)} 


2 


Il 
—2 ' ik ra ' ' 
{(G,}L “pk <G,>'{6,}|) - L “pk, ) {G.} Gey dn{e,} (eae 


1 
=) . 
- po,c ub - Gig IG)” dn{8.} + hz fic; }<G,>dn{@,} 


il IL 
-hz ,! {G,}<G,>dn{e,} + ub “(pP) E {G, }dn 


IL 


= ppic, )) (G;}<G,>dn{8] 


U2 





il 
{G.}L “p0_<G,> Wel L “pd. P {G.} ere dn{@,} (C.19) 


i 1 
-1 1 = 1 
- pul " (G,}<G,>'dn{9,} - L “ (pu) ) (6, }<G,>an {63} 


= JE 


- R 20 4 


1 it 
0, f (G,}<G,>dn{6 5} = P {6 }<G,>dn{6,} 


0 


The first term in each of the above expressions is a boundary 
term which permits incorporation of natural boundary condi- 
tions. Implementation of the boundary conditions is presented 
mimoection C.l.b. The coefficients in equations C.17, C.18, 
and C.19 are comprised of variable dependent properties, and 
were taken as the average value of the properties over an 
element. In the limit, as the elements get smaller (i.e., 
n+ -“), the average values of the coefficients converge to 
the exact values. 

Inspection of expressions C.17, C.18 and C.19 shows 


the four operators, 


f (G,}'<G,>'dn (Gaza) 
{ (6, }<G,>"dn (c.21)) 
i (6; }<G,>dn (C227) 
f {G,}dn (eee) 


Io 





To formulate these operators, the global shape function, 


G,, was defined on the local level by, 


a =) (i) 
Gy 2 ©) (a2) 


where ory and G5 were defined by, 


(l - =) for — in element (e) 
e 
oa = (G22'5) 
0 for — not in element (e) 
= for &— in element (e) 
2 
a.) (C.26) 


0 for — not in element (e) 


and de is the length of the eth element. The (+) notation 


in expression C.24 means that G. is the union of a and 
. The local shape functions (i.e., the elements) have 
the following properties, 
QR 
. Sj) om) a 
(i) ) Ge" Gy = Q ee ee a (e227) 
ie. = 3 
ae (e) _ k _ 
(ri) g. (ns) = Dae = (E223) 
O if i # j 


154 





Having defined the local shape functions, the ele- 
mental matrix operators contributing to the global matrix 


Operators C.20 through C.23 are, 


; eal 
f{G,}'<G.>tan + = (C.29) 
0 J e 
=i 1 
; ay 1 
f (Gj }<G,>' ~ * (C. 30) 
Se 
, 2 1 
f{G,}<G.> + = (C.31) 
g iL) 22 
A ; 
fied + + (6532) 
: 1 


The derivations of these operators are presented in Section 
C.l.d of this appendix. 
b. Implementation of Boundary Conditions 

Having formulated the system matrices for the field 
equations, treatment of the boundary conditions will now be 
discussed. Each field equation is considered individually. 

Peer ocrous Solid fransfer Equation 

The third set of porous solid heat transfer 

boundary conditions (i.e., expressions III.57 and III.58) will 
only be considered since the first and second sets are sub- 


sets of the third. The third set of boundary conditions for 


15> 





the porous solid are, 


36 ees 
J ll 4 “4 io 
(k +k) ne = hy (T, i> + ge (TQ ay n = QO (C355) 
*(k +k = erent = cetsers) n= 1 (C. 34) 
r 2° C5 CS 6x n ; 
Since the first term in expression C.1l7 is, 
-~_ t 
{G,}L “(1-p) (k+k_) <G,> ro) (C325) 
or in analogous form, 
2 (ep) (kK +k.) (C. 36) 
e rx on : 


Natural boundary conditions, C.33 and C.34, may be directly 
substituted into equation C.36. The response dependent 
parameters, hie hos and Toe changing with time, are evaluated 
continuously within the integration routine. Thus, the boun- 
dary conditions are incorporated in the system matrices as 
follows. 


(ly) =, alee he added to the stiffness matrix A(t) 


ig 
at location Ay 1 


C2) 2, = Galery ee T,,-oe (Te-T.) |: added to the excitation 
vector F(t) at location Fy 
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added to the stiffness matrix Ve, 
ate Location A 
3n-2,3n-2 


(3) -L*(1-p)h,: 


(4) L*(1-p) (hj T,-oe (Te-T2) |: added to the excitation 
oF F(t) ae location 


S iney, 


2. Air Heat Transfer Equation 
The first and second set of boundary conditions 
for the air heat transfer equation are, 


oT 
Lk — = pc, ulTi-T,) n 


i 
co) 


(exo) 


eS 0 n=l (C.38) 


Since these are both natural boundary conditions, they are 
substututed for the first term of expression C.18 at n = 0 
and n = l, respectively. The time dependent properties and 
parameters in the coefficient are evaluated continuously 
within the integration routine. The boundary conditions are 


implemented by adding, 


tL 


(1) +3 DA Ca teal : < the stiffness matrix A(t) at location 
SO Ae - 
ez) ppc uL TT : to the excitation vector F(t) at 
e lOecat Tonner 


Z 


The third set of air heat transfer boundary conditions are, 


Zh ay n= 0 (C. 39) 


ILS) 





a — — 
77 = 0 nN 1 (C.40) 
The essential boundary condition at n = 0 is imposed in the 


Galerkin equation as follows. The A, ; Tow of the A(t) 


~~ 


matrix, the Boi row and the Bi 2 column of the B(t) matrix, 
and the F., location of the excitation vector, F(t), are all 
set equal to zero. The Bs 2 location of the OE matrix 
is then set equal to one. 

3. Oxygen Transport Equation 


For the oxygen diffusion equation, the first and 


second set of boundary conditions are, 
bo, 22 = ule - 6.) n = 0 (C.41) 
9d = 9 n=l (C.42) 


Since these are natural boundary conditions, they were sub- 
stituted for the first term in expression C.19. The proper- 
ties are evaluated continuously within the integration routine. 


The boundary conditions are implemented by adding, 


(1) =puL-: to the stiffness matrix A(t) at location 
A a 
373 
(2) pub 79 : to the excitation vector F(t) at location 
- zi 
3 


ies! 





The third set of boundary conditions for the oxygen diffusion 


equation are, 


. a n = 0 (e743) 

0p = 

a 0 n=l (Gis44) 
The essential boundary condition at n = 0 is imposed in the 
Galerkin equation as follows. The A. , row of the A(t) 


matrix, the B. i row and the B. 3 column of the B(t) matrix, 
r] r] ot 


~s 


and the FP. location of the excitation vector, F(t), are all 
set equal to zero. The Ba 3 location of the ene matrix is 
then set equal to one. 

4. Surface Recession Problem 

As discussed previously, the oxygen diffusion 
equation is eliminated during the surface recession phase of 
a problem. By doing this, a reordering of the system matrices 
must take place, including the locations for applying the 
boundary conditions. The boundary conditions of the particle 
and air temperature equations are implemented the same as above 
except the indices identifying the matrix locations of the 
form, 3n-1i1, are changed to 2n-(i-l). 
c. Treatment of the Reaction Rate Term 

An exponential reaction rate term appears in both the 
porous solid heat transfer and the oxygen diffusion equation. 
The modified Gear integration routine requires calculating or 


approximating the Jacobian matrix from the system matrices, 


ao 





A(t) and B(t), In order to include the reaction rate term 

En the Jacobian matrix, thus improving the efficiency of the 
integration routine, the exponential terms were treated as 
follows. The reaction rate terms were multiplied and divided 


by the oxygen concentration (only the heat generation term 


is shown, oxygen consumption is treated identically), 


) 
G.9 (C.45) 
i i=l 7 33 


can 
3 


Rez @ 

g 
Applying the operation of expression C.14, and moving the 
reaction rate term divided by the oxygen concentration out- 
Side the integral, the Galerkin operator for the heat 


generation becomes, 


-l 
See Aaa 
g 2 0 


f {6,}<6,>ant95) (C.46) 
The coefficient of expression C.46 is treated as a time 
dependent parameter. The Galerkin operators are distributed 
into the stiffness matrix at locations, 3n-2,3n for heat 
generation, and 3n,3n for oxygen consumption. During the 
surface recession phase, reaction rate is treated ina 
different manner as discussed in Section III.G. 
d. Derivation of the FEM Operators 
In the section on the finite element formulation, 


the following four differential operators were identified, 
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1 


f (Gi}<Gi>dn (oa 
0 
1 
f (G5 }<Gs>dn (C.48) 
1 
i] (G5 }<G,>dn (C.49) 
0 
1 
f {G,}dn (C250) 
0 


where the G. are the global basis functions. These operators 
are constructed on the element level by introducing the 
corresponding element basis functions, g;- The global and 
element basis functions are related by, 

(1-1 


G, = gf ae ae (eo) 


where oat and g5 are defined by 
(l - =) for € in element (e) 
e 


sna = (C.52) 


0 for &—€ not in element (e) 


vee 





e for € in element (e) 
(ee) _ e 
35 = (Gros) 


0 for — not in element (e) 


and aes is the length of the (e)th element. 

The derivation of the local elemental matrices 
according to the Galerkin method for the global operations 
proceeds as follows: 


Bor @perator C.47, 


global local 
of 
il le 
f {Gj }<Gi>dn ul f <g} g4>d&]  (C.54) 
0 0 
chp 


Noting that, 


ces eee 


es 
g3 = (Ceo) 
= 


the elemental matrix becomes 


1 Ra J 
; (=) =) 1-1 
e = e 1 
f dg - a (C257) 
0 2 
ae 2 ee 
-(=) (=) — wt 
= e 


say 





For operator C.48, 


global Local 
g 
1k we 4 
) (6, }<G}>dn ut f <gy 94>dE] (C.58) 
32 


Substituting in the local shape functions gives, 


i i 
j <(- +) (=) >dé (C.59) 
= e 


(cae 2S 11 
e e e e &, I 
i Sa ey (C.60) 
0 2 2 
- (25) (25) -1 1 
22 22 
e 
For operator C.49, 
global local 
g 
iL Oe 7 
) recs — <Qy G5>dé] (C.61) 
32 


6S 





Substituting in for the local shape functions 


E 
e e ¢ ¢ 
“ <(1 = 5) (>) >dé (C.62) 
E 2 e 
(==) 
e 
the elemental matrix becomes 
2 2 
eee Sere 1 
Q (l-2 7+ =>) Ee = on ik > 
e =) Xe, e a oe 
f ; Eo ae (C.63) 
0 2 
ae an eee 1 
Sm aye Uh SB ean ay) qe 
e e 2 
e e 
mer Operator C.50, 
global local 
i 2 2 
f {G,}dn On, ag] (C.64) 
0 
0 |e, 


Substituting in the local shape function, and integrating, 


the expression becomes 


f dé + = (C.65) 
S 1 
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This last operator is used for the excitation vector as 


described in the FEM formulation. 


2. SHOOTING METHOD 

The solution of equation III.12 is obtained by the 
shooting method. A general discussion of the shooting method 
with examples is presented by Gerald [55]. The method is 
based on iterative solutions of equation III.12 until a solu- 
tion is reached which also satisfies the boundary conditions, 


P (0) Along with the boundary condition 


P., and P(1l) = Po 
at n = 0, an initial estimate to dP/dn is specified, and 
equation III.12 is integrated from n = 0 to n = 1 using 


Euler's method. An approximate Newton's method is then used 


as follows, 


dp; _ a . 
ait. Shy ee 
ee) 8 eas SOCS 88 
ae i-l ae i-2 


to provide a better approximation of dP/dn. The procedure is 
repeated until the solution has converged. A solution satis- 
fying equation III.12 and its associated boundary conditions 

is calculated at each time step. The current values of the 
properties are used and dp /ot is approximated by linear 
interpolation using the current and past values of O° The 

do /ot term was neglected for the first two integrations of a 
problem. The reason for this is associated with the difficulty 


of specifying a reasonable set of initial conditions for the 


Js 





temperature and oxygen diffusion equations. When equation 
III.12 was treated as an initial condition problem, the 
pressure gradient in a small region of the porous medium would 
approach zero for some initial conditions of temperature and 
oxygen concentration (the starting heat generation rate for 
those initial conditions was large at time, t = 0). The 
pressure gradient approaching zero resulted in numerical 


instabilities for the temperature and oxygen diffusion equations. 
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